text:
# Impact of the Wiggler Coherent Synchrotron Radiation Impedance on the Beam Instability1footnote 11footnote 1Work supported by the Department of Energy contract DE-AC03-76SF00515

## I Introduction

When an electron bunch passes the curved trajectory inside a dipole, the coherent synchrotron radiation (CSR) induces an impedance Mur97 ; Derb95 ; Warnock90 ; KYNg90 . The CSR impedance is known to impact single pass bunch compressors where the beam currents are extremely high BraunPRL ; BraunPRST ; H01 ; L02  but it is also possible that CSR might cause a microwave-like beam instability in storage rings. A theory of such an instability in a storage ring has been recently proposed in Ref. SH02  where the impedance is generated by the synchrotron radiation of the beam in the storage ring bending magnets.

A similar instability may arise due to coherent synchrotron radiation in long wigglers. Many storage rings have used large wiggler systems to reduce the damping times and control the beam emittance. In particular, the high luminosity colliding beam factories DAPHNE DAPHNE , PEP-II PEPII , KEK-B KEKB , CESR-C CESRC  all have wiggler systems that can radiate significantly more power than is radiated in the arc bending magnets. Similarly, damping ring designs for future linear colliders rely on large wiggler systems to obtain the required damping rates where, again, more power can be radiated in the wiggler systems than in the arc bending magnets.

In a previous study, the wakefield and impedance in a wiggler with large parameter $`K`$ has been obtained WRS02  based on the work of Saldin et al. Saldin98 . In this paper, we study the impact of the wiggler synchrotron radiation impedance on the beam longitudinal dynamics in rings with dipoles and wigglers. We focus our attention on the limit of a large wiggler parameter $`K`$ because this is the most interesting case for many applications.

The paper is organized as follows. In Sec. II, we will briefly review the theory of the beam instability in a ring due to the CSR impedance SH02 . Using the synchrotron radiation impedance of a wiggler WRS02 , we study the beam instability problem for a ring including both dipoles and wigglers in Sec. III. In Sec. IV, we then focus on a low-frequency microwave-like CSR instability near the pipe cut-off frequency which is the primary limitation. Finally, in Sec. V, we discuss effect of the CSR instability on the design optimization of damping rings for future linear colliders. Finally, we give some discussions and conclusion in Sec. VI.

## II Longitudinal beam instability

We consider the longitudinal dynamics of a thin coasting beam. The beam can be described with a longitudinal distribution function $`\rho (\nu ,s,z)`$. The positive direction for the internal coordinate $`s`$ is pointing to the direction of the motion. The relative energy offset of a particle having energy $`E`$ with respect to the nominal energy $`E_0`$ is expressed as $`\nu =(EE_0)/E_0`$. The position of the reference particle in the beam line is $`z=ct`$ with $`c`$ equal to the speed of light in vacuum.

The longitudinal beam instability due to the synchrotron radiation impedance in a ring has been studied theoretically by Stupakov and Heifets SH02 . In this section, we briefly review and reproduce the equations we need for our study in this paper. The reader is advised to refer to reference SH02  for the details of the theory.

For a matched beam, one can use a 1-D Vlasov equation to describe the evolution of the longitudinal distribution function $`\rho (\nu ,s,z)`$.

$`0`$ $`=`$ $`{\displaystyle \frac{\rho }{z}}\eta \nu {\displaystyle \frac{\rho }{s}}`$ (1)
$``$ $`{\displaystyle \frac{r_0}{\gamma }}{\displaystyle \frac{\rho }{\nu }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\nu ^{}{\displaystyle _{\mathrm{}}^s}𝑑s^{}w(ss^{})\rho (\nu ^{},s^{},z),`$
where, the slippage factor is defined as

$$\eta =\alpha \frac{1}{\gamma ^2},$$
(2)
with $`\alpha `$ equal to the momentum compaction factor. For ultrarelativistic beams, we have $`\eta \alpha `$. Based on our definition of $`s`$, the sign convention for $`\alpha `$ is then the following. If a particle with a higher energy than the nominal energy, i.e., $`\nu >0`$, would go to the head of the bunch, then the beam line provides a negative $`\alpha `$. According to this convention, then in a simple bending magnet, we will have $`\alpha >0`$, while in a wiggler, we have $`\alpha <0`$.

The distribution function $`\rho `$ is written as a sum of the equilibrium distribution function and a perturbation, $`\rho =\rho _0(\nu )+\rho _1(\nu ,s,z)`$, with $`\rho _1\rho _0`$. We look for a perturbation as

$$\rho _1=\widehat{\rho }_1e^{i\omega z/c+iks},$$
(3)
where $`k`$ is the wavenumber, $`\omega `$ the frequency. Linearizing the 1-D Vlasov Eq. (1) leads to the dispersion relation SH02 :

$$1=\frac{ir_0cZ(k)}{\gamma }_{\mathrm{}}^{\mathrm{}}𝑑\nu \frac{d\rho _0/d\nu }{\omega +ck\eta \nu },$$
(4)
which determines the existence of a solution as in Eq. (3). Here, $`r_02.82\times 10^{15}`$ m is the classical electron radius and $`\gamma `$ is the Lorentz factor. In addition,

$$Z(k)=_0^{\mathrm{}}𝑑sw(s)e^{iks},$$
(5)
is the CSR impedance and the wake Green function $`w(s)`$ describes the interaction of two particles due to the synchrotron radiation where $`w(s)0`$ for $`s>0`$ while $`w(s)=0`$ for $`s<0`$. The positive values of $`w(s)`$ correspond to the energy loss and the negative values imply the energy gain.

We assume that the initial energy distribution function is a Gaussian, i.e., $`\rho _0=n_0/(\sqrt{2\pi }\nu _0)\times \mathrm{exp}(\nu ^2/2\nu _0^2)`$, where $`n_0`$ is the linear density, i.e., the number of particles per unit length. In this case, Eq. (4) can be rewritten as,

$$1=\frac{iZ(k)\mathrm{\Lambda }}{\sqrt{2\pi }k}_{\mathrm{}}^{\mathrm{}}𝑑p\frac{pe^{\frac{p^2}{2}}}{\mathrm{\Omega }\pm p}$$
(6)
where, $`\mathrm{\Lambda }=n_0r_0/(|\eta |\gamma \nu _0^2)`$; $`\mathrm{\Omega }=\omega /(ck|\eta |\nu _0)`$ and $`p=\nu /\nu _0`$. The upper (lower) sign in Eq. (6) refers to the case of a positive (negative) $`\eta `$. We refer the reader to Ref. SH02  for further details.

## III Storage Rings with Wigglers

As mentioned, a number of storage rings utilize long damping wigglers where a significant fraction of the energy loss per turn is emitted. As concrete examples, we study the NLC main damping ring Andy02 , the TESLA damping ring Deck9901 , and the KEK ATF prototype damping ring ATF02 . Parameters are given in Table 1.

In our paper, we use the steady state CSR impedance and assume a distributed impedance model. For a dipole, the steady state CSR impedance is Mur97 ; Derb95

$$Z_D(k)=iA\frac{k^{1/3}}{R^{2/3}},$$
(7)
with

$$A=3^{1/3}\mathrm{\Gamma }\left(\frac{2}{3}\right)(\sqrt{3}i1).$$
(8)
The wiggler impedance $`Z_W(k)`$ is computed in Ref. WRS02 . Hence, the total impedance is then

$$Z(k)=Z_D(k)\frac{\mathrm{\Theta }R}{C}+Z_W(k)\frac{L_W}{C},$$
(9)
where, $`R`$, $`\mathrm{\Theta }`$, $`L_W`$ and $`C`$ are the dipole bending radius, the total bending angle, the wiggler total length and the damping ring circumference given in Table 1 respectively.

To study the instability, we numerically solve the dispersion relation Eq. (6). Since we are assuming a coasting beam model, we only consider instability wavelengths short compared to the bunch length.

In Fig. 1, the imaginary part of $`\mathrm{\Omega }`$ is plotted as a function of the instability wavenumber for the NLC main damping ring using the parameters listed in Table 1. At low frequencies, the dipole CSR impedance dominates while at shorter wavelengths the wiggler CSR impedance is usually more important. In the region where Im($`\mathrm{\Omega }`$)$`<0`$, the beam is stable and this is true for all regions except at the longest wavelengths as shown in the inset. This low frequency instability will be discussed in the following section.

Similar calculations were made for the TESLA “dog-bone” damping ring, where the total bending angle is about $`\mathrm{\Theta }=10\pi /3`$, and the KEK ATF prototype damping ring. For the design current in the TESLA damping ring, the impedance from the dipoles and wigglers will not drive an instability while the results for the ATF are very similar to those of the NLC damping ring shown in Fig. 1.

One interesting effect can be seen at the wiggler radiation fundamental frequency and the odd harmonics where the ring is actually stabilized by the wiggler impedance. This is opposite to the single-pass behaviour through a wiggler where there is an instability at the wiggler fundamental frequency usually refered to as the FEL instability. The effect will be discussed further in Appendix A but is a direct result of the positive sign of the momentum compaction in the ring.

To illustrate this, Fig. 2 is a plot for the NLC damping ring with identical wiggler and arc parameters but the momentum compaction is assumed to be opposite in sign. Here, one can see that at low frequencies the growth is slightly lower. In contrast, the wiggler impedance makes the system less stable at the wiggler fundamental unlike the case with positive momentum compaction. Furthermore, if the magnitude of the momentum compaction is reduced, the system will become unstable similar to that illustrated in Appendix B. This is similar to the FEL instability however, it is also noted in Appendix A that our theory does not correctly treat the full FEL instability.

## IV Low-frequency CSR instability

As seen in Figs. 1 and 2, the instability is most important at relatively low frequency. The longest instability wavelength that is possible is determined by two effects: first, the instability wavelength has to be short compared to the bunch length and, second, the vacuum chamber causes an exponential suppression of the synchrotron radiation at wavelengths $`\lambda `$ greater than the ‘shielding cutoff’ Warnock90

$$\lambda _c4\sqrt{2}b(b/R)^{1/2}.$$
(10)
Here, $`R`$ is the dipole bending radius, and $`b`$ is the vacuum chamber half height. The numerical coefficient of Eq. (10) assumes that the vacuum chamber is made up of two infinitely wide plates. Different cross sections give different numerical factors SK02 . Given the previous discussion, the threshold will be the lowest at the smaller of the bunch length or the ‘shielding cutoff’ wavelength.

For the NLC main damping ring, we find that perturbations with wavelengths $`\lambda >3.5`$ mm are not stable due to the dipole CSR impedance alone. Adding the CSR impedance from the wiggler causes perturbations with wavelengths $`\lambda >2.6`$ mm to be unstable. In Fig. 3, the threshold particle number is plotted as a function of the perturbation wavelength. Next, in Fig. 4, the growth rate, which is defined as the inverse of the time needed for the perturbation to grow by a factor of $`e`$

$$\frac{1}{\tau }\mathrm{Im}(\mathrm{\Omega })ck|\eta |\nu _0.$$
(11)
is plotted versus the perturbation wavelength. Based on the parameters in Table 1 and Eq. (10), the ‘shielding cutoff’ wavelength is computed to be $`\lambda _c4.9`$ mm. At this cutoff wavelength, the threshold currents and growth time are summarized in Table 1. The growth time is significantly faster than the synchrotron period, in agreement with the analysis for a microwave-like instability.

It is clearly seen that the threshold current decreases as we approach the longer wavelength perturbations, however, the growth rate is not monotonic. This is the result of two opposite mechanisms: one is the energy modulation growth due to the CSR impedance, and the other is the Landau damping. For very long wavelength perturbation, Landau damping effect is small, hence, we could expand the denominator of the integrand in Eq. (4) to get

$$\omega =\sqrt{ic^2n_br_0\eta kZ(k)/\gamma },$$
(12)
which indicates that $`\tau ^1`$ Im($`\omega `$) $`\sqrt{kZ(k)}`$, i.e., the growth rate will increase for a shorter wavelength perturbation, since $`Z(k)k^\epsilon `$, with $`\epsilon >0`$. This is shown in Fig. 4: when we approach from the long wavelength perturbation to the short wavelength perturbation, the growth rate increases.

However, this process stops when Landau damping becomes effective. The finite energy spread in the beam will produce a phase mixing due to the slippage and this will destroy density modulation due to the CSR induced energy modulation. The Landau damping due to the phase mixing is more serious for short wavelength perturbations. This damping can be seen in the second term of Eq. (1), or more clearly, in the denominator of the dispersion relation in Eq. (4). This is demonstrated in Fig. 4, the growth rate finally decreases when we approach very short wavelength, and eventually, the system becomes stable.

For the KEK ATF prototype damping ring, the cut-off wavelength would be about $`\lambda _c3.1`$ mm according to Eq. (10). In Fig. 5, the threshold particle number is plotted as a function of the perturbation wavelength. Taking the dipole CSR impedance alone, for the single bunch charge, the instability sets in for perturbations with wavelengths $`\lambda >2.8`$ mm. Adding the wiggler CSR impedance, the electron beam would be unstable for perturbations with wavelengths $`\lambda >1.9`$ mm. Other results are summarized in Table 1.

It is interesting to note that in both the NLC and the ATF damping rings roughly twice as much synchrotron radiation power is emitted in the wiggler as in the arc dipoles. However, the instability threshold is not dramatically impacted by the additional radiation in the wiggler and decreases by less than a factor of two in each case. This arises because of the very different low frequency behavior of the impedances.

We can use a simple scaling analysis to understand this. According to Eq. (7), the dipole CSR impedance scales as $`Z_D(k)k^{1/3}`$. Similarly, the low frequency behavior of the wiggler impedance is WRS02

$$Z_W(k)=\pi k_w\frac{k}{k_0}\left[1\frac{2i}{\pi }\mathrm{log}\left(\frac{k}{k_0}\right)\right]$$
(13)
which is accurate enough for practical calculation in $`k[0,0.1k_0]`$, where $`k_0`$ is the wiggler fundamental radiation wavenumber:

$$k_0=2\gamma ^2k_w/(1+K^2/2)$$
(14)
and the wiggler parameter $`K`$ is approximately $`K93.4B_w\lambda _w`$, with $`B_w`$ the peak magnetic field of the wiggler in units of Tesla and $`\lambda _w`$ the period in meters. Thus, the wiggler CSR impedance scales as Re($`Z_W(k)`$)$`k`$, and Im($`Z_W(k)`$)$`k\mathrm{log}(k)`$, which have a weaker scaling than $`Z_D(k)k^{1/3}`$ when $`k0`$, i.e., the CSR impedance from the wiggler is weaker than that from the dipole when we approach low frequency region.

Hence, it suggests an optimization of the damping ring design where a larger fraction of radiation is emitted in the wigglers as will be discussed subsequently.

## V Damping Ring Optimization

Damping rings are used in linear colliders to attain the very small transverse emittances that required to attain the desired small spot sizes at the collision point. In the damping rings, the emission of synchrotron radiation damps the phase space volume of the injected beams towards the equilibrium values. An $`e`$-folding reduction of the transverse phase space occurs after the emission of synchrotron radiation energy comparable to the beam energy.

There are usually three conflicting requirements on the damping rings: first, the equilibrium emittances must be small which usually requires weak bends and strong focusing; second, the dynamic aperture must be large to accept the injected beams and this is difficult to attain in a strong focusing storage ring; and, third, the damping rates must be fast so that the incoming phase space and transients can be damped before the beam are extracted and accelerated to the collision point.

To achieve these conflicting requirements most damping ring designs have been forced to include long damping wigglers where a significant fraction of the energy loss per turn is emitted. Of course, the ring characteristics can be greatly modified by changing the ratio of the damping due to the wiggler versus that due to the bending magnets. However, historically the wiggler systems have been a concern because they are a significant source of non-linearity and because they generate intense radiation which may directly or indirectly impact the beam.

Other work has improved the confidence in the dynamic aperture predictions and has shown that properly designed wigglers should not impose a limitation marco03 . Thus, we can further optimize the ring designs based on quantitative evaluations of the effect of the wiggler radiation fields.

To study the ring optimization, we use expressions derived in Ref. ER02  for the damping ring parameters. The primary constraint is the damping times. Because all linear collider designs presently only consider flat beams with $`\sigma _x\sigma _y`$, the tightest requirement is on the vertical damping time which, using Eq. (11) of Ref. ER02 , can be written

$$\tau _y\frac{(2.89\times 10^{12}\mathrm{kG})C}{|B_a|\gamma ^2(1+F_w)c}\frac{|\mathrm{\Theta }|}{2\pi },$$
(15)
where $`B_a`$ is the arc bending field, $`C`$ is the ring circumference, and $`\mathrm{\Theta }`$ is the total bending angle of the arcs. In addition, $`F_w`$ is a parameter equal to the ratio of the damping due to the wiggler over the damping due to the arc bending magnets which can be written:

$$F_wI_{2w}/I_{2a}.$$
(16)
Here, $`I_{2a}`$ and $`I_{2w}`$ are the second synchrotron integrals calculated over the arcs and the wigglers respectively, which can be written

$$I_{2a}=\left|\frac{\mathrm{\Theta }}{2\pi }\frac{2\pi B_a}{(B\rho )}\right|=\left|\frac{\mathrm{\Theta }}{R}\right|;I_{2w}\frac{L_wB_w^2}{2(B\rho )^2}.$$
(17)
where $`(B\rho )`$ is the standard energy dependent magnetic rigidity, $`(B\rho )0.017\gamma `$ (kG m), and $`B_a`$, $`B_w`$, and $`L_w`$ are the magnetic field of the bending dipole, the peak wiggler field, and the wiggler length.

According to Eq. (15), we can scale the bend field as

$$B_a\frac{1}{(1+F_w)}$$
(18)
for a given ring size and damping time and, given a maximum wiggler field, the required wiggler length can be found from Eqs. (17). At this point, we can calculate the other ring parameters. According to the dispersion relation in Eq. (6), we will need the momentum compaction factor $`\alpha `$, the energy spread $`\nu _0`$, and the bunch rms length $`\sigma _z`$. Using Eqs. (37), (39) and (40) from Ref. ER02 , the scaling can be written:

$`\alpha (1+F_w)^{5/3};`$ $`\nu _0\sqrt{{\displaystyle \frac{B_a+B_wF_w}{1+F_w}}};`$
$`\sigma _z`$ $``$ $`\nu _0\sqrt{\alpha }.`$ (19)
Examples of the scaling of these parameters versus $`F_w`$ can be found in Figs. 11 and 13 of Ref. ER02  for an earlier version of the NLC damping ring. Using these relations, we estimated parameters for NLC ring with $`F_w=0`$ and $`F_w=6.5`$. These parameters are summarized in Table II along with the nominal parameters where $`F_w=2.2`$.

At this point, we can calculate the impedance. Based on our previous studies, we expect the CSR impedance from the wiggler to be small compared to the CSR impedance of dipoles in the low frequency region of interest. Hence, the threshold is essentially determined by the dipole CSR impedance. According to the impedance in Eqs. (9) and (7), and the relation of $`B_a`$ and $`F_w`$ in Eq. (15), the total impedance should scale as

$$Z(k)Z_D(k)RR^{1/3}(1+F_w)^{1/3}.$$
(20)
at a fixed wavenumber $`k`$. Clearly, this is an increasing function with $`F_w`$ however, according to the dispersion relation in Eq. (6), the threshold will scale as:

$$N_t\frac{\alpha \nu _0^2\sigma _z}{Z(k)}(1+F_w)^{13/6},$$
(21)
where it was assumed that $`B_aB_w`$, so $`\nu _0`$ is roughly constant. Thus, it is expected that increasing $`F_w`$ will increase the threshold significantly.

In addition, the cutoff wavelength $`\lambda _cR^{1/2}B_a^{1/2}(1+F_w)^{1/2}`$. Hence, as $`F_w`$ increases and $`B_a`$ decreases, the cutoff wavelength will decrease. Since, the threshold is lowest at long wavelength, the threshold will increase as $`F_w`$ is increased due to the shift in the cutoff wavelength as well as the change in the ring parameters.

The above scaling analysis has been compared to numerical calculations for the NLC damping ring. The results are given in Fig. 6 where we plot the threshold particle number as a function of the CSR instability wavelength for $`F_w=0`$ (dashes), $`2.2`$ (solid) and $`6.5`$ (dots). In addition, the cutoff wavelength for the three cases is indicated with vertical lines using the same notation. Clearly, the threshold increases rapidly as $`F_w`$ is increased as expected. The approximately quadratic scaling given above is followed as $`F_w`$ is increased from the nominal value of 2.2 to 6.5 but the dependence is weaker as $`F_w`$ is decreased to 0 because of the change in the beam energy spread which was assumed constant. Based on these results a new version of the NLC damping ring has been designed with an $`F_w6.5`$ which has a significantly higher threshold for both this CSR instability as well as the conventional microwave instability driven by the vacuum chamber impedance newNLCDR .

## VI Discussion and conclusion

In this paper, we studied a possible beam instability due to coherent synchrotron radiation from dipole and wiggler magnets in a storage ring. This instability may be important in many high luminosity factories such as DAPHNE, PEP-II, KEK-B, and CESR-C as well as the damping rings of future linear colliders. We used the results to analyze the stability of the damping ring designs and found that a possible instability due to the CSR will exist in the NLC damping ring design and the ATF prototype damping ring at KEK. We then discuss optimization of the rings and the pros and cons of increasing the damping fraction due to the wigglers. This led to a redesign for the NLC damping which increases the threshold to more than four times the nominal charge per bunch.

As we discussed above, in our calculation, we used steady state impedance for both the dipoles and wigglers. In the future, it may be important to consider transient effects in our calculation. In addition, the shielding cutoff of the CSR should be calculated as part of the impedance and not treated as an ad hoc limit. This is also important because we assumed the same shielding cutoff around the ring but in practice it will differ in the bending magnets and the wigglers. Finally, the bunch length is the same order of the instability wavelength, while in the theory, we use a coasting beam model. In the future, the finite bunch length should be taken into consideration. All of these topics are under study and will be reported in a separate paper.

###### Acknowledgements.

The authors thank A.W. Chao, P. Emma, S.A. Heifets, M. Ross, M. Woodley of Stanford Linear Accelerator Center, S. Krinsky, J.B. Murphy, J.M. Wang of National Synchrotron Light Source, Brookhaven National Laboratory, and K.J. Kim of Advanced Photon Source, Argonne National Laboratory for many discussions. Work was supported by the U.S. Department of Energy under contract DE-AC03-76SF00515.

## Appendix A Growth rate at the FEL frequency and the odd harmonics

In this appendix, we discuss the growth rates calculated at the wiggler fundamental frequency, frequently refered to as the FEL frequency, and at the odd harmonics. The wiggler CSR long-range wakefield is WRS02

$$G(\zeta )=\frac{1}{2\zeta }+\frac{1}{2\zeta }\underset{n=0}{\overset{\mathrm{}}{}}\left[J_n\left(\frac{2n+1}{2}\right)J_{n+1}\left(\frac{2n+1}{2}\right)\right]^2\mathrm{cos}(4(2n+1)\zeta ).$$
(22)
To the leading order in $`k`$, this gives us the CSR impedance as

$`Z(k)`$ $`=`$ $`\pi k_w{\displaystyle \frac{k}{k_0}}\left\{1{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left(\left[J_n\left({\displaystyle \frac{2n+1}{2}}\right)J_{n+1}\left({\displaystyle \frac{2n+1}{2}}\right)\right]^2H[k(2n+1)k_0]\right)\right\}2ik_w{\displaystyle \frac{k}{k_0}}\mathrm{log}\left({\displaystyle \frac{k}{k_0}}\right)`$ (23)
$`+`$ $`{\displaystyle \frac{1}{2}}ik_w{\displaystyle \frac{k}{k_0}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left\{\left[J_n\left({\displaystyle \frac{2n+1}{2}}\right)J_{n+1}\left({\displaystyle \frac{2n+1}{2}}\right)\right]^2\mathrm{log}\left({\displaystyle \frac{[k(2n+1)k_0]^2[k+(2n+1)k_0]^2}{k_0^4}}\right)\right\},`$
where, $`H(x)`$ is the Heaviside step function, i.e., $`H(x)=1`$ for $`x>0`$, $`H(0)=1/2`$, and $`H(x)=0`$ for $`x<0`$.

At the FEL frequency and the odd harmonics, the impedance has a logarithmically divergent imaginary part, Im$`(Z(k))\mathrm{log}(|\mathrm{\Delta }k|/k)`$, with $`\mathrm{\Delta }kk(2n+1)k_0`$, and the real part is finite. We write $`Z=RiL`$, where $`R>0`$ is the resistance and $`L>0`$ is the inductance, and $`LR`$. As a first estimate, we neglect the resistance and write $`Z=iL`$ only. In the dispersion relation Eq. (6), then $`|\mathrm{\Omega }|L1`$. Now because of the factor $`\mathrm{exp}(p^2/2)`$, $`\mathrm{\Omega }`$ is huge compared with $`p`$ in the $`p`$-nonvanishing integral region. Hence the integrand could be expanded. Keeping the first non-zero term, the integral is simplified as

$$_{\mathrm{}}^{\mathrm{}}𝑑p\frac{pe^{\frac{p^2}{2}}}{\mathrm{\Omega }\pm p}\sqrt{2\pi }\frac{1}{\mathrm{\Omega }^2},$$
(24)
for large $`\mathrm{\Omega }`$ to get

$$\mathrm{\Omega }^2=\pm \frac{L\mathrm{\Lambda }}{k}.$$
(25)
The upper (lower) sign in Eq. (25) refers to the case of a positive (negative) $`\eta `$. For the rings which we study in this paper, they all have $`\eta >0`$ and we take the upper sign. Hence, Im$`(\mathrm{\Omega })=0`$, i.e., the growth rate will be zero and the system is stable at the FEL frequency and all the odd harmonics. This agrees with the conventional knowledge that a ring is stable with a purely inductive impedance and a positive momentum compaction or a purely capacitive impedance and a negative momentum compaction.

## Appendix B Single wiggler

In this appendix, we study the beam instability in a single wiggler. From FEL theory, we would expect instability at the wiggler fundamental radiation frequency and all the harmonics. In addition, we need to look at a microwave-like low-frequency instability.

In a wiggler, the slippage factor is

$$\eta =\frac{\left(1+\frac{K^2}{2}\right)}{\gamma ^2}.$$
(26)
i.e., the slippage factor is negative, therefore, we choose the lower ($``$) sign in Eq. (6). In Fig. 7, we plot the imaginary part of $`\mathrm{\Omega }`$ as a function of the instability wave number for the wigglers in the NLC and TESLA damping rings Andy02 ; Deck9901 . For the region where Im($`\mathrm{\Omega }`$)$`<0`$, the beam is stable. We find that, for the nominal current, the system is stable within the entire interesting region where the perturbation wavelength $`\lambda <3`$ cm, except at the wiggler fundamental radiation wavelength and the third harmonic.

In fact, according to Eq. (26), we have $`\eta <0`$ and we need take the lower sign in Eq. (25). When we approach the fundamental and the third harmonic frequency, $`L`$ becomes large, and Eq. (25) gives us a large Im($`\mathrm{\Omega }`$), i.e., a large growth rate as we find in Fig. 7. However, according to Eq. (25), we will have Im$`(\mathrm{\Omega })I^{1/2}`$, where $`I`$ is the peak current of the electron bunch; while in the FEL theory, the growth rate scales as Im$`(\mathrm{\Omega })I^{1/3}`$ BPN84 . To correctly calculate the FEL instability from the wake impedance approach, the electron density distribution function has to be taken at the retarded time in the standard Vlasov equation SK03 . This then leads to a cubic dispersion relation same as in the FEL theory, and results in the cubic root scaling in the growth rate: Im$`(\mathrm{\Omega })I^{1/3}`$.

Nevertheless, here we further use the FEL theory to study the instability at the FEL frequency and the third harmonic. We use the simulation code TDA TW89  to calculate the instability growth length at the fundamental frequency and a newly developed code JHW02  to calculate the instability at the third harmonic. The wiggler and beam parameters used are listed in Table 1. In the NLC damping wiggler, the growth length at the FEL frequency is about 50 m assuming that the phase errors between wiggler sections is small. Fortunately, the separation between the 2-meter long wiggler sections will further reduce the gain. In addition, any modulation that is induced at the fundamental wavelength of $`10\mu `$m will be smeared out by the energy spread and the momentum compaction as the beam is transported around the remainder of the ring. A similar conclusion is found for the TESLA damping ring wigglers and the instability at the FEL frequency and the odd harmonics are not expected to be a limitation in either case.


format: text

text:
# Jet Fragmentation due to a Quark/Diquark Pick-up in High Energy Heavy Ion Collisions

## Acknowledgements

This work is partly supported by the US DOE grant. We thank B.Kopeliovich, C.M.Ko, R.J.Fries and T. Sjöstrand for helpful comments.


text:
# Stochastic Stabilization of Cosmological Photons

The stability of photon trajectories in models of the Universe that have constant spatial curvature is determined by the sign of the curvature: they are exponentially unstable if the curvature is negative and stable if it is positive or zero. We demonstrate that random fluctuations in the curvature provide an additional stabilizing mechanism. This mechanism is analogous to the one responsible for stabilizing the stochastic Kapitsa pendulum. When the mean curvature is negative it is capable of stabilizing the photon trajectories; when the mean curvature is zero or positive it determines the characteristic frequency with which neighbouring trajectories oscillate about each other. In constant negative curvature models of the Universe that have compact topology, exponential instability implies chaos (e.g. mixing) in the photon dynamics. We discuss some consequences of stochastic stabilization in this context.

One of the fundamental questions concerning the dynamical properties of photons in the Universe is whether their trajectories are stable or unstable. This strongly influences both the images of distant objects as well as fluctuations in the cosmic microwave background (CMB). In cosmological models in which the spatial geometry of the Universe has constant curvature $`K`$ the photon trajectories (geodesics) are stable if $`K>0`$, in which case neighbouring trajectories oscillate about each other with a characteristic frequency determined by $`K`$, and exponentially unstable if $`K<0`$, in which case neighbouring trajectories diverge exponentially quickly with a Liapunov exponent determined by $`|K|`$.

It is obvious, however, that the Universe is not exactly homogeneous and isotropic: matter is not uniformly distributed but is organized into galaxies, clusters of galaxies, and even superclusters of galaxies . Consequently, the spatial curvature cannot be constant, but fluctuates. Our purpose here is to investigate the influence of these fluctuations on the stability of the photon trajectories. It might be thought that random fluctuations in the curvature would be a source of instability. We show that this is not in fact the case; their influence is a stabilizing one in that the curvature is renormalized by an additional (positive) factor related to the fluctuation amplitude and length-scales. The mechanism responsible for this stochastic stabilization is analogous to the one that stabilizes the Kapitsa pendulum: a pendulum whose pivot is forced to move periodically or stochastically along the vertical .

Given the fact that for constant curvature models to be consistent with recent observations of the CMB the absolute value of the mean curvature must be small – the density of energy in the Universe is within about $`1\%`$ of its critical value  – stochastic stabilization is capable of dominating the stability balance. For example, an order of magnitude estimate suggests that it is relevant in the present cosmological epoch. The difference between a constant curvature model with $`K`$ close to zero and one in which stochastic stabilization dominates is quantifiable in terms of the stability frequency with which neighbouring trajectories oscillate about each other.

Our results apply to both open and compact geometries. There has, however, recently been considerable attention focused on negative curvature models that have compact topology . Photon dynamics in compact spaces is recurrent and hence, by virtue of the exponential instability, strongly chaotic (e.g. mixing) when the curvature is constant and negative. The fact that fluctuations in the cosmic microwave background (CMB) have a distribution close to gaussian has been related to Berry’s random wave model  for wave modes in chaotic systems  and the phenomenon of scarring in these wave modes  has been linked with anisotropic structures in the CMB and in the distribution of galaxies . We discuss some implications of stochastic stabilization in this context.

We begin by deriving the appropriate form of the geodesic deviation equation – the equation that governs the separation of neighbouring trajectories. Since the aim of this Letter is to establish the principle of stochastic stabilization for photon trajectories, the calculations we report contain the essential ingredients for the effect, but ignore many additional, comparatively weaker phenomena present in the early Universe. With this in mind, we begin with an expanding isotropic homogeneous (Friedmann) cosmology perturbed by density fluctuations with non-relativistic velocities, neglecting perturbations of vector and tensor character such as gravitational waves, and pressure fluctuations, for example caused by relativistic neutrinos. In the coordinate system called the ‘conformal Newtonian’ or ‘longitudinal’ gauge , with the spatial variables in Robertson-Walker form, such a spacetime is described by the metric

$$\begin{array}{c}ds^2=R^2(\tau )\{(1+2\mathrm{\Phi })d\tau ^2\frac{12\mathrm{\Phi }}{[1+\frac{K}{4}(x^2+y^2+z^2)]^2}(dx^2+dy^2+dz^2)\}.\hfill \end{array}$$
(1)
Here $`R(\tau )`$ is the scale factor or spatial curvature radius of the Universe, $`\tau =R^1𝑑t`$ is conformal time, $`\mathrm{\Phi }<<1`$ is the Newtonian gravitational potential, $`K`$ is the dimensionless spatial curvature (e.g. $`K=1`$ corresponds to a hyperbolic geometry), $`(x,y,z)`$ are comoving coordinates expanding at the same rate as the Universe, and units are chosen in which Newton’s constant $`G`$ and the speed of light $`c`$ are equal to unity. The paths of photons through the spacetime are null geodesics described by the equation

$$\frac{d𝐧}{d\tau }=2_{}\mathrm{\Phi }$$
(2)
where $`𝐧`$ is the photon direction, and $`_{}`$ is the gradient in comoving coordinates perpendicular to $`𝐧`$. Note that since $`\mathrm{\Phi }`$ and its derivative are small, the total change of $`𝐧`$ is small, and so the transverse derivative can be replaced by the derivative transverse to the observed (final) direction of the photon. The separation of two closely spaced photons propagates according to the geodesic deviation equation  which under our approximations leads to

$$\frac{d^2}{d\tau ^2}\left(\begin{array}{c}\xi \\ \eta \end{array}\right)=\left(\begin{array}{cc}K+2\mathrm{\Phi }_{\xi \xi }& 2\mathrm{\Phi }_{\xi \eta }\\ 2\mathrm{\Phi }_{\xi \eta }& K+2\mathrm{\Phi }_{\eta \eta }\end{array}\right)\left(\begin{array}{c}\xi \\ \eta \end{array}\right)$$
(3)
where $`\xi `$ and $`\eta `$ are the components of the separations orthogonal to $`𝐧`$, and there are contributions arising from both the spatial curvature $`K`$ and second derivatives of the gravitational potential in directions orthogonal to $`𝐧`$, namely tidal forces.

These tidal forces are, in principle, entirely characterized by a complete knowledge of the fluctuations in the matter density. We shall consider them to be random functions of position (and possibly time). As it moves under their influence, a photon will be seen to be acted on by a time-varying force which fluctuates rapidly (because the speed of light is large and the length scale of the density fluctuations is small compared to $`R`$) and randomly (i.e. stochastically) with zero mean. We thus replace the spatially random tidal force terms in (3) with a time-dependent stochastic perturbation.

In order to illustrate the qualitative behaviour of the solutions of the resulting class of equations, we consider first the analogous one-component case:

$$\frac{d^2u}{d\tau ^2}=(k+Af(\tau ))u.$$
(4)
Here $`A`$ is a control parameter and $`f(\tau )`$ is a stochastic forcing function, which we take to have zero mean. This has a mechanical analogy: it describes a pendulum with a vertically moving pivot in the limit of small oscillations. In this case $`u`$ denotes the angular displacement from the vertical and the forcing term describes the height of the pivot. This problem was studied in detail by Kapitsa . The constant gravitational force represented by $`k`$ ($`k<0`$ corresponds to an up-turned pendulum and $`k>0`$ to a down-turned one) plays the role of the smooth geometry in the case of the unperturbed cosmology, while the motion of the pivot corresponds to the metric perturbations induced by fluctuations in the matter density.

We shall be interested in the case when the frequency $`\omega `$ characterizing the fluctuations in $`f(\tau )`$ is large. Stabilization has been proved in the limit $`\omega \mathrm{}`$ using Liapunov exponent techniques (see, for example, ). The dependence of the stability on $`\omega `$ can be deduced by separating $`u`$ asymptotically into fast and slow components : $`u(\tau )=u(\tau )+u_f(\tau )`$, where $`\mathrm{}`$ denotes a local time average over scales large compared to $`\omega ^1`$. Here $`u(\tau )`$ is the slow component and, as $`\omega \mathrm{}`$, $`u_f`$ is small and varies rapidly. Substituting into (4) and averaging, we have

$$\frac{d^2u}{d\tau ^2}=kuAf(\tau )u_f(\tau ).$$
(5)
Subtracting (5) from (4) gives, as $`\omega \mathrm{}`$,

$$\frac{d^2u_f}{d\tau ^2}Af(\tau )u(\tau ).$$
(6)
This equation can be integrated directly, treating $`u(\tau )`$ as a constant, leading to

$$\frac{du_f}{d\tau }Au(\tau )_{\tau _1}^\tau f(t)𝑑tAu(\tau )v(\tau ),$$
(7)
where $`\tau _1`$ is chosen so that $`v=0`$, and

$$u_f(\tau )Au(\tau )_{\tau _2}^\tau v(t)𝑑tAu(\tau )x(\tau ),$$
(8)
$`\tau _2`$ being chosen so that $`x=0`$. Substituting into (5) and computing the integral implicit in the average by parts, we obtain

$$\frac{d^2u}{d\tau ^2}(k+A^2v^2(\tau ))u.$$
(9)
It follows from (9) that if $`A^2v^2<k`$ then $`u`$ grows exponentially as $`\tau \mathrm{}`$, but if $`A^2v^2>k`$ then $`u`$ is linearly stabilized. This effect, and the accuracy of (9), are illustrated by numerical simulations, the results of which are represented in Figures 1, 2 and 3. In these simulations we took $`k=1`$ and

$$f(\tau )=\underset{n=1}{\overset{N}{}}\mathrm{sin}(\omega _n\tau +\varphi _n)$$
(10)
with $`N=100`$, frequencies $`\omega _n`$ chosen at random uniformly from $`[120,600]`$ and phases $`\varphi _n`$ chosen at random uniformly from $`(\pi ,\pi ]`$. Thus $`v^2=\frac{1}{2}_n\frac{1}{w_n^2}\frac{1}{1440}`$ and so for $`A<\sqrt{1440}38`$ the expectation is that $`u`$ grows like $`\mathrm{exp}(\tau \sqrt{1A^2/1440})`$, while for $`A>\sqrt{1440}`$ it oscillates with frequency $`\sqrt{A^2/14401}`$.

The above analysis generalizes directly to the geodesic deviation equation (3). This may be written as the stochastic matrix equation

$$\frac{d^2}{d\tau ^2}\left(\begin{array}{c}\xi \\ \eta \end{array}\right)=\left(\begin{array}{cc}K+Af_{11}(\tau )& Af_{12}(\tau )\\ Af_{12}(\tau )& K+Af_{22}(\tau )\end{array}\right)\left(\begin{array}{c}\xi \\ \eta \end{array}\right).$$
(11)
The variables $`\xi `$ and $`\eta `$ can be separated into slow and fast components, as for $`u`$. Following the steps that led from (4) to (9) we find, in the high-frequency limit, that

$$\frac{d^2}{d\tau ^2}\left(\begin{array}{c}\xi \\ \eta \end{array}\right)(K\text{I}+A^2\text{v}^2)\left(\begin{array}{c}\xi \\ \eta \end{array}\right)$$
(12)
where I is the $`2\times 2`$ identity matrix and v is the matrix with elements $`v_{ij}=^\tau f_{ij}(t)𝑑t`$, such that $`v_{ij}=0`$. The stability of the trajectories thus depends on the eigenvalues of the matrix $`K\text{I}+A^2\text{v}^2`$; specifically, assuming $`v_{11}=v_{22}`$, they are unstable if $`A^2(v_{11}\pm v_{12})^2<K`$, and if $`A^2(v_{11}\pm v_{12})^2>K`$ they are stable. Note that $`v_{11}`$ and $`v_{12}`$ are of the order of $`\omega ^1`$, and so the balance between the terms in the stability criterion depends on the size of $`A^2/\omega ^2|K|`$. In particular, since $`\omega `$ is considered large here, if $`K<0`$ the perturbations must be large enough to cause the curvature to be positive in places for stabilization to occur. This is consistent with the fact that the geodesics on compact surfaces where the curvature is everywhere negative (but not necessarily constant) are strongly chaotic.

Let us make an order of magnitude estimate to determine the cosmological epoch when the stochastic term is observable. Since the curvature is small, it will always be dominated by an observable stochastic term (see below). Consider that the Universe consists of a volume fraction $`(R_{}/R)^3`$ of clusters of galaxies of size $`L1`$Mpc and mass $`M10^4`$Mpc in relativistic units. The number of clusters within a ball of radius $`R`$ is thus $`(R_{}/L)^3`$. The scale factor $`R`$ is set to the scale of the visible Universe at present, roughly $`10^4`$Mpc. $`R_{}`$ is the distance scale at which voids between the galaxies began to form, at a redshift of order 10 (i.e. $`R_{}10^3`$Mpc). The mean time (in terms of $`\tau `$) for a photon to move between clusters is of the order of $`LR^2/R_{}^3`$. This gives $`\omega R_{}^3/(LR^2)`$, which is large up to the present, but will fall to less than one in a few times the current age of the Universe, as the galaxies move apart. Note that $`\omega `$ is assumed large in the calculations described above.

The tidal forces $`\mathrm{\Phi }_{\xi \xi }`$ are, up to factors of $`4\pi `$, of the order of the dimensionless density $`MR^2/L^3`$, from Poisson’s equation when a photon passes near a cluster; $`Av`$ is typically this quantity integrated over a mean free time, $`MR^4/L^2R_{}^3`$. $`|K|0.01`$ in these units. The stochastic stabilization will be observable if $`Av1`$, and will dominate the curvature if $`Av/|K|1`$ (a weaker condition). Taking this together with the estimate for $`\omega `$, we conclude that stochastic stabilization is potentially observable when $`(L^2/MR_{})^{1/4}R/R_{}(R_{}/L)^{1/2}`$; that is, for a period starting soon after the galaxies formed and due to end a few times the present age of the Universe from now. The magnitude of each deflection is determined by $`M/L10^4`$, so this effect should be observable on angular scales of arcseconds or smaller . From this perspective it should be easier to observe in distant quasars, which are almost pointlike, than the almost homogeneous cosmic background radiation.

The stochastic stabilization of photon trajectories has several interesting consequences for cosmology. In particular when $`K0`$ it is observable in the characteristic frequency with which neighbouring orbits oscillate about each other, corresponding to a periodic focusing and defocusing of distant images.

Given the recent focus on compact models of the Universe , we devote some concluding remarks to this case (although, as already noted, the mechanism is independent of topology). The simplest models of the Universe that go beyond having spatially constant negative curvature incorporate static fluctuations. Such models are to a large extent still unrealistic, because the curvature is unlikely to be static over the timescales associated with photon recurrence. Nevertheless, they illustrate most clearly the implications of stabilization, some of which we now list.

As a result of stabilization, the photon dynamics will generically possess both regular (stable) and irregular (chaotic) components, rather than being fully chaotic. While wave modes in fully chaotic systems are believed to have a gaussian value distribution, it is well established in the context of quantum chaos that regular trajectories lead to a quantifiable non-gaussian component whose precise form depends on the size and position of the regular regions in phase space . This will in turn give rise to a non-gaussian component in the CMB, via the connection described in . There is, in addition, likely to be a second non-gaussian component due to the photon trajectories close to bifurcation (i.e. close to making the transition between instability and stability). It was established in  that bifurcations give rise to non-gaussian fluctuations in wave modes, quantified by characteristic scaling exponents in the wavelength dependence of the moments of their value distribution in the short-wavelength limit. In situations when all of the generic bifurcations contribute, the moment exponents take on universal values, which were calculated in . These non-gaussian statistics are also likely to feed through to the fluctuations in the CMB in the models under consideration here.

Bifurcations have a strong influence on the scarring of wave modes as well. In fully chaotic systems wave modes may be scarred by periodic orbits . This effect is dramatically enhanced when the orbits undergo bifurcation , giving rise to superscars. It has been suggested by Levin and Barrow, in the context of constant negative curvature models, that scars may give rise to anisotropic structures in the CMB . This mechanism would therefore be strongly enhanced by stabilization. Levin and Barrow  have also put forward the idea that the anisotropic structures associated with scars may show up in the distribution of matter. Regions of regular motion and superscars would significantly amplify this effect. In particular, stable islands typically have a fractal distribution, and this may provide an explanation, within the context of their theory, for the fractal hierarchy of structures seen in galaxies and clusters of galaxies. Stable orbits also strongly inhibit the rate of mixing, giving rise to intermittency, and so may play a role in the suggested links between mixing and preinflationary homogenization .

We gratefully acknowledge very helpful discussions with Dr Gopal Basak and Professors Sir Michael Berry and Mark Birkinshaw.


text:
# Chapter 1 Numerical Likelihood Analysis of Cosmic Ray Anisotropies

## Chapter 1 Numerical Likelihood Analysis of Cosmic Ray Anisotropies

### Abstract

A numerical likelihood approach to the determination of cosmic ray anisotropies is presented which offers many advantages over other approaches. It allows a wide range of statistically meaningful hypotheses to be compared even when full sky coverage is unavailable, can be readily extended in order to include measurement errors, and makes maximum unbiased use of all available information.

### 1. Introduction

The search for anisotropies in the cosmic ray distribution (and of other objects – see for example ) is of increasing interest with the advent of the Pierre Auger Observatory (for a review of the physics up to this time, see ). The issues involved are, however, subtle and complicated by limited statistics at the highest energies and nonuniform sky coverage. The first full-sky search for UHECR anisotropies using a standard approach is presented in this conference . This paper presents a maximum likelihood approach which is well-adapted to further studies in anticipation of much larger data sets.

### 2. Expansion in Spherical Harmonics

It is convenient to define a set of real spherical harmonics $`\{\psi _{\mathrm{},m}\}`$ for $`\mathrm{}=0,1,2,\mathrm{}`$ and $`m=\mathrm{},(\mathrm{}1),\mathrm{},0,\mathrm{},(\mathrm{}1),\mathrm{}`$ by

$$\psi _{\mathrm{},m}(\theta ,\varphi )=\{\begin{array}{cc}k_{\mathrm{}}^{|m|}P_{\mathrm{}}^{|m|}(\mathrm{cos}\theta )\mathrm{cos}(m\varphi )\hfill & \text{for }m=\mathrm{},\mathrm{},1,\hfill \\ k_{\mathrm{}}^0P_n^0(\mathrm{cos}\theta )\hfill & \text{for }m=0\text{,}\hfill \\ k_{\mathrm{}}^mP_{\mathrm{}}^m(\mathrm{cos}\theta )\mathrm{sin}(m\varphi )\hfill & \text{for }m=1,\mathrm{},\mathrm{}\hfill \end{array}$$
(1)
which are orthonormal with respect to the usual measure $`\mathrm{sin}(\theta )d\theta d\varphi `$ and integration over $`\theta `$ from 0 to $`\pi `$ and $`\varphi `$ from 0 to $`2\pi `$, the $`P_{\mathrm{}}^m`$ are associated Legendre polynomials of the first kind, and the normalization constants $`k_m^{\mathrm{}}`$ are:

$$k_0^0=\frac{1}{\sqrt{2\pi }},k_{\mathrm{}}^0=\sqrt{\frac{2\mathrm{}+1}{4\pi }},\mathrm{and}k_{\mathrm{}}^m=\sqrt{\frac{2\mathrm{}+1}{2\pi }\frac{(\mathrm{}m)!}{(\mathrm{}+m)!}}$$
(2)
The natural measure of anisotropy for a spherical distribution is in terms of these spherical harmonics, as each $`\mathrm{}`$ labels, in a coordinate-independent fashion, just how much of each irreducible $`SO(3)`$ representation is present.

In a perfect world with infinite statistics and complete sky coverage there are now many possible approaches to estimating how much of each of these components is present in a distribution, or, better, what is the likelihood that a given function $`f(\theta ,\varphi )`$ with Fourier-Legendre expansion $`f(\theta ,\varphi )=_{\mathrm{}=0}^{\mathrm{}}_{m=\mathrm{}}^{\mathrm{}}a_{\mathrm{},m}\psi _{\mathrm{},m}`$ representing the probability density of sources gives rise to the observed distribution $`g(\theta ,\varphi )`$. The coefficients can be extracted from the usual integral

$$a_{\mathrm{},m}=_0^{2\pi }_0^\pi f(\theta ,\varphi )\psi _{\mathrm{},m}(\theta ,\varphi )\mathrm{sin}\theta d\theta d\varphi $$
(3)
Of course this gives no measure of what sort of error should be associated with the determined values of each coefficient. Alternative approaches are to fit for the coefficients by minimizing some $`\chi ^2`$-like quantity like

$$_0^{2\pi }_0^\pi \frac{\left|g(\theta ,\varphi )_{\mathrm{}=0}^{\mathrm{}}_{m=\mathrm{}}^{\mathrm{}}a_{\mathrm{},m}\psi _{\mathrm{},m}(\theta ,\varphi )\right|^2}{\sigma ^2}\mathrm{sin}\theta d\theta d\varphi $$
(4)
with $`\sigma `$ a suitable measure of error, or to compute and maximize a corresponding likelihood that the hypothesized distribution parametrized by the $`a_{\mathrm{},m}`$ gives rise to the the observed distribution $`g`$. Of course this is statistically unreasonable for a finite number of sources (i.e. to provide an infinite number of coefficients!), In general a decision must be made to truncate the expansion at some value of $`\mathrm{}`$ to make the sum finite, but now the general phenomenon of aliasing risks that, for example, a fit allowing for $`\mathrm{}=0,1`$ might give misleading results for observed data which is drawn from a purely $`\mathrm{}=2`$ distribution, say.

A more serious problem is that should part of the sky be unobserved, there is now no way to calculate anything! This is not a trivial point. An attempt to find anisotropies based only on observations in the Northern hemisphere with zeroes inserted for the whole Southern hemisphere would be wildly in error if some simple extrapolation were made to the unobserved region of the sky – especially if there were something bright and as-yet undetected in the South! The real challenge is to say something statistically meaningful with the data that is actually available. Clearly, observing only the Northern hemisphere and seeing a good degree of isotropy should increase one’s net belief in overall isotropy of the full sky, while leaving open the possibility of a staggeringly bright or empty sky in the South. One approach is to try to make functions which would be orthonormal over the observed region of the sky, but it’s not clear that this has much physical meaning as it elevates lack of acceptance to a status comparable to the $`SO(3)`$ invariance of space embodied in the $`\psi _{\mathrm{},m}`$. The following is a proposal for what seems to make good statistical sense.

### 3. A Likelihood Proposal

Based on the above observations, we make the following proposal: keep the spherical harmonics as always with the (necessarily) truncated Fourier-Legendre expansion and construct an unbinned likelihood  function in which one clearly specifies which values of $`\mathrm{}`$ are included in the sum. An unbinned likelihood function, as described below, automatically makes maximum use of all detected information, allows for measurement errors to be included easily, and is easy to implement numerically. Most importantly, however, the likelihood is not to be normalized as it stands. Rather, we take the Bayesian approach to likelihood which says that likelihoods give us ways to update our prior experimental data or guesses in light of new information. This will mean that we are able to present results on various hypotheses about data without bias, and with the easy inclusion of other data from the same, or other experiments.

To be concrete, we specialize here to the case where $`g`$ is a sum of delta functions representing sources $`i`$ of unit intensity at $`\theta _i,\varphi _i`$ and later discuss how one can treat the case of these being at uncertain locations, or taking into account other properties such as intensity, energy, composition, etc.. These will appear as natural generalizations to the approach described.

The (unnormalized!) likelihood $`L(\{\mathrm{}\}|\theta _i,\varphi _i)`$ that the measured $`\theta _i,\varphi _i`$ arise from $`f(\theta ,\varphi )=_{\mathrm{}}^{\mathrm{}}_{m=\mathrm{}}^{\mathrm{}}a_{\mathrm{},m}\psi _{\mathrm{},m}`$ where the sum over $`\mathrm{}`$ is specified according to whatever hypothesis is being tested (i.e. just taking $`\{\mathrm{}\}=\{0,1\}`$ allows for uniform and dipole contributions and no others, while $`\{\mathrm{}\}=\{2\}`$ would be pure quadrupole) is

$$L(\{\mathrm{}\}|\theta _i,\varphi _i)=\frac{f(\theta _i,\varphi _i)}{_{\mathrm{}}_{m=\mathrm{}}^{\mathrm{}}da_{\mathrm{},m}f(\theta ,\varphi )}$$
(5)
where the integral in the denominator is over all the parameters that can vary and over the range in which the parameters are allowed to vary. An important caveat is that for $`f`$ to represent a sensible probability distribution it should never be negative, and this must be checked. Two approaches are possible: one is to restrict the domain over which coefficients range so that the function is strictly positive (this is not actually very difficult in practice since distributions are often nearly uniform with small fluctuations superimposed) or to take as a probability distribution function some positive function of $`f`$ in place of $`f`$ above. In astronomy  it is not uncommon to see $`\mathrm{exp}(f)`$. The choice ultimately represents the unavoidable presence of some (often hidden) assumption about what a sensible prior (i.e. in the absence of data) is for the likelihood – a point to which we return later.

By construction then $`L(\{\mathrm{}\}|\theta _i,\varphi _i)`$ is normalized so that its integral over all the parameters that can vary (the coefficients $`a_{\mathrm{},m}`$ included in the truncated Fourier-Legende expansion, and the total likelihood $`L_{TOT}`$, which is a function of those $`a_{\mathrm{},m}`$ is

$$L_{TOT}(\{\mathrm{}\})=\underset{i}{}L(\{\mathrm{}\}|\theta _i,\varphi _i)$$
(6)
This is a relative likelihood, and while no absolute normalization is possible (nor should it be!) if part of the sky is unobserved, it is now very useful for two types of calculation. If one has a prior expectation for the distribution of the $`a_{\mathrm{},m}`$ (which might be that they are all a priori equally likely) then $`L_{TOT}`$ can be multiplied by this and the product treated as a normalized likelihood distribution for the $`a_{\mathrm{},m}`$ themselves. This is, of course, potentially dangerous, but does allow one to see how the new data (the measured $`(\theta _i,\varphi _i)`$) should cause one to revise earlier beliefs. Such a likelihood can be maximized with respect to the $`a_{\mathrm{},m}`$ and if not enough data is present to test the hypothesis (for example, more coefficents to fit for than data points) the fit will respond by just not converging (that is to say, there will be flat directions in the likelihood as a function of the parameters meaning that one can’t decide) - the beauty of this approach is that it is, by construction, correct. When results are obtained, one automatically gets the values of the parameters, and their entire likelihood distributions from which errors (which need not even be Gaussian) can be extracted. One can even do things like fix some parameters to those given by a favoured theory and then repeat the fitting and then obtain likelihoods for the correctness of that theory.

More objectively, one can compute relative likelihoods in which the prior drops out, so that it is reasonable to ask (even in the absence of full sky coverage!) what the relative likelihood $`L_{REL}`$ is of pure dipole distribution to one admitting uniform, dipole and quadrupole components: $`L_{REL}=\frac{L_{TOT}(\{1\})}{L_{TOT}(\{0,1,2\})}`$ If one wants statements made to be about all energies over $`x`$, then one just uses the data points with energies over $`x`$. Similarly, data can be selected by composition and (even relative) anisotropies be searched for as functions of composition, energy, time, etc. Extensions to uncertainty in direction are trivial to include: simply divide a given event into a large number, $`M`$, say of subevents distributed appropriately and count each in the likelihood with a weight $`1/M`$ – this numerically convolves this uncertainty with the likelihood, and the size of any additional errors introduced by the procedure can be studied by varying $`M`$.

### 4. References

1. Anchordoqui L., Paul T., Reucroft S., Swain J., 2002 arXiv:hep-ph/0206072.

2. Anchordoqui, L., Hojvat, C., McCauley, T. P.,Paul, T.C., Reucroft, S., Swain, J. D., Widom, A., these proceedings and arXiv:astro-ph/0305158.

3. Edwards, A. W. F. 1972. Likelihood. Cambridge University Press, Cambridge.

4. Evans, N. W., Ferrer, F., Sarkar, S., Astropart. Phys. 17 (2002) 319.

5. Peebles, P. J. E., Ap. J. 185 (1973) 413; Hauser, M.G., Peebles P. J. E., Ap. J. 185 (1973) 757.


text:
# 1 Introduction

## 1 Introduction

In 2013, the author asked the following quesiton titled “Is there a ‘mathematical’ definition of ‘simplify’?” on the popular mathematics website MathOverflow.net :

“Every mathematician knows what ‘simplify’ means, at least intuitively. Otherwise, he or she wouldn’t have made it through high school algebra, where one learns to ‘simplify’ expressions like $`x(y+x)+x^2(y+1+x)+3(x+3)`$. But is there an accepted rigorous ‘mathematical’ definition of ‘simplify’ not just for algebraic expressions but for general expressions, which could involve anything, like transcendental functions or recursive functions? If not, then why? I would think that computer algebra uses this idea.”

The answers there indicated that even though every mathematician knows intuitively what “simplify” means, there is still no universally accepted definition of “simplify”. In fact, one of the answers (by Henry Cohn) indicated that “In full generality, there provably isn’t any method for complete simplification”. (He was referring to elementary functions of a real variable.) In this paper, we shall give a simple and plausible definition of “simplify” in terms of the computational complexity of integer functions. We shall also use this definition to show that there is no deterministic and exact algorithm which can compute the permanent of an $`n\times n`$ matrix in $`o(2^n)`$ time.

## 2 A definition of “simplify”

Consider the following definition of “simplify”:

Definition: An algebraic expression (recursive or non-recursive) for a function $`f:`$ cannot be simplified if there is no other algebraic expression for $`f`$ which can be computed faster.

For example, the expression $`xw+yz+xz+yw`$ can be simplified to $`(x+y)(w+z)`$, since computing $`(x+y)(w+z)`$ takes only one multiplication and two additions, while computing $`xw+yz+xz+yw`$ takes four multiplications and three additions. And we can also see clearly that the expression $`(x+y)(w+z)`$ cannot be simplified.

As another example, let $`f:`$ be the function which satisfies the recursive formula, $`f(n)=f(n1)+1`$ and $`f(0)=0`$. This recursive formula can be simplified to $`f(n)=n`$, since computing the recursive formula for $`f`$ takes $`\mathrm{\Theta }(n)`$ time, while computing the formula $`f(n)=n`$ is trivial. And the formula $`f(n)=n`$ clearly cannot be simplified.

And let $`f:`$ be the function which satisfies the recursive formula, $`f(n)=f(n1)+f(n2)`$ and $`f(1)=f(2)=1`$, the Fibonacci sequence. This recursive formula can be simplified, since it is possible to prove that $`f(n)`$ equals $`\varphi ^n/\sqrt{5}`$ rounded to the nearest integer, where $`\varphi =(1+\sqrt{5})/2`$, which can be computed exponentially faster than the recursive formula can be computed .

## 3 Computing the permanent of a matrix

Let $`A=(a_{ij})`$ be a matrix of integers. The permanent of $`A`$ is defined as:

$$\text{perm}(A)=\underset{\sigma S_n}{}\underset{i=1}{\overset{n}{}}a_{i\sigma (i)},$$

where $`S_n`$ is the symmetric group . The fastest known deterministic and exact algorithm which computes the permanent of a matrix was first published in 1963 and has a running-time of $`\mathrm{\Theta }^{}(2^n)`$ . It is still considered an open problem by the mathematics and computer science community whether this time can be beaten. Now consider the following theorem and proof, which we shall discuss afterwards:

Theorem: There is no deterministic and exact algorithm which can compute the permanent of an $`n\times n`$ matrix in $`o(2^n)`$ time.

Proof: For any row $`i`$, the permanent of matrix $`A`$ satisfies the recursive formula

$$\text{perm}(A)=\underset{j=1}{\overset{n}{}}a_{ij}\text{perm}(A_{ij}^\mathrm{\#})$$

and $`\text{perm}([a_{11}])=a_{11}`$, where $`A_{ij}^\mathrm{\#}`$ is the $`(n1)\times (n1)`$ matrix that results from removing the $`i`$-th row and the $`j`$-th column from $`A`$. This formula cannot be simplified, so the fastest algorithm for computing the permanent of a matrix is to apply this recursive formula to matrix $`A`$. Since this involves recursively evaluating the permanent of $`\mathrm{\Theta }(2^n)`$ submatrices of $`A`$, each corresponding to a subset of the $`n`$ columns of $`A`$, we obtain a lower bound of $`\mathrm{\Theta }(2^n)`$ for the worst-case running-time of any deterministic and exact algorithm that computes the permanent of a matrix.∎

At first, this proof makes sense intuitively, but if one thinks about it a little more, one might become skeptical, since one could argue the same for the determinant of a matrix, that there is no deterministic and exact algorithm which can compute the determinant of an $`n\times n`$ matrix in $`o(2^n)`$ time (which is known to be false) - for any row $`i`$, the determinant satisfies the recursive formula

$$\text{det}(A)=\underset{j=1}{\overset{n}{}}(1)^{i+j}a_{ij}\text{det}(A_{ij}^\mathrm{\#})$$

and $`\text{det}([a_{11}])=a_{11}`$, which is almost the same as the recursive formula for the permanent of a matrix.

However, there is a big difference between the two recursive formulas: There are negative signs in the formula for the determinant, so it is not inconceivable that one might be able to cancel most of its terms out, if one is clever. And in fact this is the reason why it is possible to compute the determinant of a matrix in polynomial-time: If one performs elementary row operations on matrix $`A`$ with pivot $`a_{11}0`$, converting it to a matrix $`B`$ with zeroes in the last $`n1`$ entries of column 1, then the determinant of $`A`$ will equal the determinant of $`B`$ and we will also obtain a simpler formula for the determinant:

$$\text{det}(A)=a_{11}\text{det}(B_{11}^\mathrm{\#}).$$

This trick ultimately leads to a polynomial-time algorithm for computing the determinant of a matrix, if one applies it recursively to the matrix $`B_{11}^\mathrm{\#}`$, exchanging rows when necessary.

However, in the case of the permanent of a matrix, no trick like this is possible, since there are only positive signs in its formula. To gain some insight as to why this is so, consider the following analogy: Suppose we want to subtract two large positive numbers with a tiny difference, say $`a=12,345,678,907`$ and $`b=12,345,678,903`$. One could compute $`a`$ minus $`b`$ by applying the normal subtraction procedure that one learns in elementary school to each digit of these two numbers, but one does not have to do this; if we let $`c=12,345,678,900`$, then we will obtain the same answer by computing $`(ac)`$ minus $`(bc)`$, which amounts to subtracting only the last digits of each number, $`7`$ minus $`3`$. But there are no short-cuts like this for adding $`a`$ and $`b`$, since none of their digits can be cancelled out. And for this same reason, it is possible to cancel out lots of terms in the formula for the determinant but not in the formula for the permanent, as the elementary row operations which are performed on matrix $`A`$ when computing its determinant via the algorithm described above are analogous to subtracting $`c`$ from both $`a`$ and $`b`$.

But then one might ask, “The proof above said ‘This formula cannot be simplified’. But how can I be sure of this?” The answer to this question is that we know that the above recursive formula for the permanent cannot be simplified, because we have tried every possible way to simplify it and saw that each way fails: To be specific, we tried to multiply the factors, $`a_{ij}`$ and $`\text{perm}(A_{ij}^\mathrm{\#})`$, of the summands together, but we failed since the two factors are completely independent from one another. And we tried adding the summands together, but we also failed since the factors $`a_{ij}`$ found in each summand are completely independent from one another and are also completely independent from each $`\text{perm}(A_{ij}^\mathrm{\#})`$; furthermore, we found that since $`\text{perm}(A_{ij}^\mathrm{\#})`$ is different in each term, it is impossible to use the distributive law to decrease the computational complexity of the recursive expression. And finally, we noticed that the row choice of $`i`$ is irrelevant in the recursive formula for the permanent, so no choice of $`i`$ is better than any other choice. What other things are there to try that could possibly make the expression simpler? Nothing, since we have already considered every mathematical operation in the recursive formula for the permanent. Therefore, the recursive formula for the permanent cannot be simplified, i.e., it has the best computational complexity of any algebraic expression for the permanent of a matrix.

This type of reasoning is not new or foreign; it is essentially the same type of reasoning that a high school math student uses to simplify algebraic expressions. Also note that only if one is careful in one’s analysis and considers every possible way to simplify an algebraic expression can one prove that an algebraic expression indeed cannot be simplified; merely claiming that an algebraic expression cannot be simplified does not make it so. But sometimes it is so obvious that an algebraic expression cannot be simplified that writing down a full explanation of this is unnecessary. Also, it turns out that one can use similar reasoning to prove that there is no deterministic and exact algorithm which solves the Traveling Salesman Problem in polynomial-time .

## 4 Conclusion

While everyone in the mathematics community understands intuitively what “simplify” in mathematics means, there is still no universal definition of “simplify”. In this paper, we have defined “simplify” in terms of the computational complexity of an integer function and have shown that this definition can be used to prove that there is no deterministic and exact algorithm which can compute the permanent of an $`n\times n`$ matrix in $`o(2^n)`$ time.


text:
# Cosmic Microwave Background and Supernova Constraints on Quintessence: Concordance Regions and Target Models

## Abstract

We perform a detailed comparison of the Wilkinson Microwave Anisotropy Probe (WMAP) measurements of the cosmic microwave background (CMB) temperature and polarization anisotropy with the predictions of quintessence cosmological models of dark energy. We consider a wide range of quintessence models, including: a constant equation-of-state; a simply-parametrized, time-evolving equation-of-state; a class of models of early quintessence; scalar fields with an inverse-power law potential. We also provide a joint fit to the CBI and ACBAR CMB data, and the type 1a supernovae. Using these select constraints we identify viable, target models for further analysis.

The precision measurement of the cosmic microwave background (CMB) by the Wilkinson Microwave Anisotropy Probe (WMAP) satellite MAPmission ; MAPresults  represents a milestone in experimental cosmology. Designed for precision measurement of the CMB anisotropy on angular scales from the full sky down to several arc minutes, this ongoing mission has already provided a sharp record of the conditions in the Universe from the epoch of last scattering to the present. In light of this powerful data Hinshaw:2003ex ; Kogut:2003et ; Page:2003fa ; Spergel:2003cb ; Verde:2003ey , we must consider anew our cosmological theories.

We aim to use the WMAP results to test cosmological theories of the accelerating Universe — to seek clues to the nature of the dark energy. Despite the absence of a direct dark-energy interaction with our baryonic world, the CMB photons provide a probe of the presence of the dark energy, complementary to the type 1a supernovae. Via the integrated Sachs-Wolfe effect on large angular scales, the geometric optics of the last-scattering sound horizon on degree scales, and the pattern of acoustic oscillations on smaller angular scales, we expect the CMB to reveal information about the dark energy density, equation-of-state, and behavior of fluctuations. These effects are illustrated in Figure 1.

In this article we test the quintessence hypothesis, that a dynamical, time-evolving, negative pressure, inhomogeneous form of energy dominates the cosmic energy density and is responsible for the cosmic acceleration Ratra:1987rm ; Peebles:1987ek ; Wetterich:fm ; Wetterich:bg ; Coble:1996te ; Caldwell:1997ii ; Turner:1998ex . To be precise, we carry out an extensive analysis of the CMB anisotropy and mass fluctuation spectra for a wide range of quintessence models. These models are: (Q1) models with a constant equation-of-state, $`w`$, including $`w<1`$; (Q2) models with a simply-parametrized, time-evolving $`w`$; (Q3) early quintessence models, with a non-negligible energy density during the recombination era; (Q4) trackers described by a scalar field evolving under an inverse-power law potential.

The suite of parameters describing the cosmological models are split into spacetime plus “matter sector” variables, $`\theta _M`$, and separate quintessence parameters, $`\theta _Q`$. The spacetime and “matter sector” of the quintessence models are specified by the parameter set $`\theta _M=\{\mathrm{\Omega }_bh^2,\mathrm{\Omega }_{cdm}h^2,h,n_s,A_S,\tau _r\}`$. In order, these are the baryon density, cold dark matter density, hubble parameter, scalar perturbation spectral index, scalar perturbation amplitude, and optical depth. In this investigation we restrict our attention to spatially-flat, cold dark matter models with a primordial spectrum of nearly scale-invariant density perturbations generated by inflation.

The quintessence parameters vary from model to model. For the simplest family of models, with a constant equation-of-state, we need only to specify $`\theta _Q=\{w\}`$. For models which feature more realistic time-evolution of the quintessence, which may include a non-negligible fraction of quintessence at early times, more parameters are required, e.g. $`\theta _Q=\{w,dw/da,\mathrm{}\}`$, to characterize the impact on the cosmology in general and the CMB in particular.

Our analysis method is as follows: (i) compute the CMB and fluctuation power spectra for a given cosmological model; (ii) compute the relative likelihood of the model with respect to the experimental data; (iii) assemble the likelihood function in parameter space to determine the range of viable quintessence models. For step (i) we use both a version of CMBfast CMBfast  modified for quintessence, as well as the newly-available CMBEASY CMBEASY . For step (ii) we supplement the WMAP data with the complementary ACBAR ACBAR  and CBI-MOSAIC data CBI\_mosaic ; CBI\_deep  (using the same bins as Ref. Spergel:2003cb ; Verde:2003ey ), in addition to the current type 1a SNe data Schmidt:1998ys ; Riess:1998cb ; Garnavich:1998th ; Perlmutter:1998np . Certain other constraints, such as the HST Key Project measurement of $`H_0`$ Freedman:2000cf  or the limit from Big Bang Nucleosynthesis on $`\mathrm{\Omega }_bh^2`$ Burles:2000zk  through the deuterium abundance measurement are satisfied as cross-checks. It is remarkable that such agreement can be found between such diverse phenomena. Our focus in the following investigation, however, is primarily on the CMB and SNe.

Q1: We have analyzed the cosmological constraints on the simplest model of quintessence, characterized by a constant equation-of-state, $`w`$. We have used the equivalence between a scalar field $`\phi `$ with potential $`V(\phi )`$ and the equation-of-state $`w`$ in order to self-consistently evaluate the quintessence fluctuations. For the range $`w<1`$ we employ a k essence model, keeping the sound speed (actually, this is $`d\omega ^2/dk^2`$) fixed at $`c_s^2=1`$. Since this model introduces only one additional parameter beyond the basic set of spacetime plus matter sector variables, we adopt a simplistic grid-based search for viable models. The acceptance criteria for the Q1 models is based on a $`\mathrm{\Delta }\chi ^2`$-test. The results of our survey of Q1 models are shown in Figures 2, 3. We have exploited the degeneracy of the CMB anisotropy pattern among models with the same apparent angular size of the last scattering horizon Huey:1998se . Hence, there is a family of models with $`\mathrm{\Omega }_bh^2=0.023`$, $`\mathrm{\Omega }_{cdm}h^2=0.126`$, $`n_s=0.97`$, and characterized by pairs $`\{w,h\}`$ having (nearly) indistinguishable CMB anisotropy patterns. The pairs $`\{w,h\}`$ are shown in Figure 2, and all represent quintessence models with $`\chi ^2=1429`$ for the WMAP temperature-temperature and temperature-polarization data — a one-parameter family of best-fit models. From this starting point, we explored over $`6\times 10^4`$ models distributed on a grid filling a 6-cylinder around the best-fit line, varying $`\{\mathrm{\Omega }_bh^2,\mathrm{\Omega }_{cdm}h^2,h,n_s,\tau _r\}`$ at intervals in $`w`$. For each model we evaluate the likelihood relative to WMAP, as well as the complementary ACBAR and CBI-MOSAIC data. The $`2\sigma `$ boundary, based on a $`\mathrm{\Delta }\chi ^2`$ test for six degrees of freedom is shown in Figure 3. We have also evaluated the constraint in the $`w\mathrm{\Omega }_m`$ plane for the combined Hi-Z/SCP type 1a supernova data set, showing the $`2\sigma `$ region based on a $`\mathrm{\Delta }\chi ^2`$ test for two degrees of freedom. Our basic conclusion from the overlapping constraint regions is that there exist concordant models with $`1.25w0.8`$ and $`0.25\mathrm{\Omega }_m0.4`$. We have identified four sample models in Table 1 for further analysis.

Our search of the parameter space for the remaining models is based on a Bayesian approach, using a Monte-Carlo Markov-Chain (MCMC) search algorithm to identify the best cosmological models. The end product is a realization of the posterior probability distribution function on the parameter space Christensen:2001gj . Our approach is similar to the procedure described in Verde:2003ey , whereby the MCMC makes a “smart” walk through the parameter space, accepting or rejecting sampled points based on a running criteria. For each of Q2-4, after some experimentation we found it practical to use four independent chains in order to monitor convergence and mixing according to the criteria of Ref. Gelman::1992 . Each such chain explored $`3\times 10^4`$ models. We found it advantageous to use a dynamical stepsize factor $`r[0.1,10]`$. In each step, all parameters are varied according to a Gaussian with width $`r\sigma _i`$, where $`\sigma _i`$ is an initial stepsize for the random walk in the direction of the $`i`$-th parameter. The stepsize factor $`r`$ is adjusted such that the chain neither rejects nor accepts too frequently. Certain parameters, such as $`w_0`$ and $`\mathrm{\Omega }_Q^{ls}`$ in the early quintessence models, have physical boundaries and it is therefore important that our MCMC does not lead to a bias in the sampled distribution close to the boundaries. As the likelihood in the MCMC approach is given by the frequency with which the chain passes through a given region of parameter space, we would be in danger of depleting the region close to the boundary if we were to simply discard any steps which cross the boundary. That is, the likelihood in the parameter region close to the boundary would be artificially suppressed. If such a (virtual) jump in a parameter direction $`i`$ occurs, we use the following strategy. We try making a gaussian-weighted step but only towards the boundary. If this fails after ten attempts, then we start over again by making a gaussian step in any direction. By trying ten times, however, it is very likely that we do make the step successfully, if we are at least one $`\sigma _i`$ away from the boundary. If, however, we are very close to the boundary already and our step would have carried us deep into the unphysical region, then it is sensible to make the step away from the current point, for the chain didn’t intend to stay close to the boundary anyhow. To give us confidence in the reliability of this algorithm, we have checked by synthetic distributions that this strategy is very well suited for sampling the likelihood, even in cases where a physical boundary cuts off the parameter space.

Q2: We have examined quintessence models with an equation-of-state that evolves monotonically with the scale factor, as $`w(a)=w_0+(1a)w_1`$. For this case, the parameters are simply $`\theta _Q=\{w_0,w_1\}`$. This parametrization has been shown to be versatile in describing the late-time quintessence evolution for a wide class of scalar field models Linder:2002et . Based on the degeneracy of models found for Q1, we expect to find a two-dimensional family of equivalent best-fit models with the same apparent angular size of the last scattering horizon, occupying a plane in the $`\{w_0,w_1,h\}`$ space. There are three ways in which this plane is pared down: Firstly we confine $`w1`$ at all times. Secondly, the transition from $`w_0`$ to $`w_0+w_1`$ takes place at low redshift, $`z1`$, so that high redshift supernovae restrict $`w_0,w_1`$ for these models. Thirdly, this parametrization allows for models in which the dark energy is non-negligible at early times, which influences the small-scale fluctuation spectrum. The first two considerations yield $`w_0<0.75`$ at the $`2\sigma `$ level, marginalizing over the suppressed five dimensional parameter space, as illustrated in Figure 4(a). There, the shapes of the contours indicate that current data can only distinguish between fast ($`dw/da0.5`$) and slow evolution of $`w(a)`$, and offer only a weak bound on $`w_1`$. However, in terms of $`\mathrm{\Omega }_Q^{ls}`$, our third consideration gives a tight upper bound on the quintessence density during recombination. As shown in Figure 4(d), $`\mathrm{\Omega }_Q^{ls}<0.03`$ at the $`2\sigma `$ level. Three target models, with significantly different equation-of-state evolution $`dw/da`$, are given in Table 2 for future investigations.

Q3: We have examined models of leaping kinetic quintessence, a scalar field evolving under an exponential potential with a non-canonical kinetic term that undergoes a sharp transition at late times, leading to the current accelerated expansion Hebecker:2000zb . At early times the field closely tracks the cosmological background with $`w=0`$ during matter domination, appearing as early quintessence Caldwell:2003vp  before undergoing a steep transition towards a strongly negative equation-of-state by the present day. The steepness of the transition in $`w`$ for a leaping kinetic model is directly connected to the equation-of-state $`w_0`$ today. Such models can therefore be characterized by the parameters $`\theta _Q=\{\mathrm{\Omega }_Q^{ls},w_0\}`$, where $`\mathrm{\Omega }_Q^{ls}`$ is the quintessence density during recombination. (A more general parameterization, allowing for independent $`w_0`$ and steepness of transition can be found in Corasaniti:2002vg .) Since $`\mathrm{\Omega }_Q^{ls}`$ is not tied as closely to the expansion rate sampled by the supernovae, compared to case Q2, the result is the weaker constraint $`\mathrm{\Omega }_Q^{ls}0.1`$, as shown in Figure 4(b). Although the limit of a cosmological constant can be approached in this model, the presence of early quintessence will then require a sharp transition in the equation-of-state in order to reach $`w1`$. In addition to the fact that such models lose the early tracking behavior and instead require some degree of fine-tuning, there is the practical consideration that the sharp transition leads to some numerical instability in our code. To avoid this problem, we restrict $`w>0.97`$ as can be seen in Figure 4(b). We believe that we have succeeded at implementing an MCMC algorithm that does not artificially distort the probability distribution near the parameter-space boundaries, based on our experimentation with synthetic distributions with boundaries. Next, because early quintessence suppresses the growth of fluctuations on small scales compared to large scales, we find that comparable fluctuation spectra can be achieved by making a trade-off between $`n_s`$ and $`\mathrm{\Omega }_Q^{ls}`$. As shown in Figure 4(e), slight suppression of small-scale power can be accomplished either by a tilt towards the red, $`n_s<1`$, or a rise in $`\mathrm{\Omega }_Q^{ls}`$. Since the effect of early quintessence on the small-scale fluctuation power spectrum closely mimics a running spectral index, we have not introduced $`\mathrm{d}n_s/\mathrm{d}\mathrm{ln}k`$ as an additional parameter, which would be highly degenerate with $`\mathrm{\Omega }_Q^{ls}`$ Caldwell:2003vp . We expect that improved measurements of the second and third acoustic peaks will tighten the constraint on $`\mathrm{\Omega }_Q^{ls}`$ and sharpen the degeneracy in the $`n_s\mathrm{\Omega }_Q^{ls}`$ plane. Target models, with significantly different values of early quintessence abundance, are given in Table 3.

Q4: Finally, we have examined tracker models of quintessence. Inverse-power law (IPL) models are the archetype quintessence models with tracking property and acceleration Ratra:1987rm ; Zlatev:1998tr ; Steinhardt:nw . The potential is given by $`V\phi ^\alpha `$, where the constant of proportionality is determined by $`\mathrm{\Omega }_Q`$. In certain supersymmetric QCD realizations of the IPL Masiero:1999sq , $`\alpha `$ is related to the numbers of color and flavors, and can take on a continuous range of values $`\alpha >0`$. For $`\alpha 0`$, however, inverse-power law models behave more and more like a cosmological constant. Using earlier data, $`\alpha `$ has been constrained to $`\alpha <1.4`$ Balbi:2000kj ; Baccigalupi:2001aa , although $`h=0.65`$ has been fixed in those analyses. Keeping $`h`$ free, a more conservative value of $`\alpha <2`$ Doran:2001ty  was inferred. From our analysis, we see that the $`2\sigma `$ bound has not changed dramatically. The $`2\sigma `$-bound with $`\alpha 12`$ is consistent with values of $`h`$ within the range determined by the HST, as seen in Figure 4(c). In Figure 4(f) we plot the likelihood contours in the $`\mathrm{\Omega }_mh^2\alpha `$ plane: our results agree with the best fit at $`\mathrm{\Omega }_mh^2=0.149`$ for $`\alpha =0`$ or $`w=1`$, but show a tolerance for a wider range for $`0\alpha 2`$. That is, the additional degree of freedom in $`\alpha `$ means that the matter density for the IPL model is not as well-determined from the peak position Page:2003fa  as compared to the $`\mathrm{\Lambda }`$ model. However, to maintain the peak at $`\mathrm{}=220`$, we observe that $`\mathrm{\Omega }_mh^2`$ decreases slightly as $`\alpha `$ increases. We might have inferred the results for the IPL based on the constant equation-of-state models: pairs of $`\{\alpha ,h\}`$ can equivalently determine a family of models with degenerate CMB anisotropy patterns, since the differences in the late-ISW for this model compared to Q1 make only a small contribution to the overall $`\chi ^2`$. Furthermore, since IPL quintessence can be modeled by the appropriate choice of the Q2 parameters, then improved sensitivity to $`dw/da`$ is required to tighten the constraints here. While $`\alpha `$ is tightly constrained, our take-away is that IPL models with $`0.25\mathrm{\Omega }_m0.4`$ remain viable. Target models are given in Table 4.

This work provides a capsule summary of the viable quintessence dark energy models, based on two of the tightest constraint methods, using the CMB and SNe. Our study advances beyond past investigations Baccigalupi:2001aa ; Brax:2000yb ; Corasaniti:2001mf ; Bean:2001xy ; Hannestad:2002ur ; Bassett:2002fe ; Jimenez:2003iv ; Barreiro:2003ua  by treating a wide class of quintessence models with the powerful weight of the WMAP data. We have considered four classes of models which cover the most basic quintessence scenarios, including a versatile parametrization, as well as the best motivated and most realistic scenarios based on our current understanding of particle physics. Absent from our survey are k essence models, and dark energy models with a coupling to other matter fields or gravity. In the former case, since the sound speed of fluctuations varies with time in these models, however, we can make a simple distinction with quintessence models with an underlying scalar field, wherein the propagation speed is equal to the speed of light. (See Refs. Erickson:2001bq ; DeDeo:2003te  for analysis of these models with respect to CMB anisotropy.) For the latter case we refer to Ref. Amendola:2003eq  for specific coupled models. We also note that the mass fluctuation power spectrum is an important cosmological constraint which we have omitted at this stage, primarily because it constrains energy density rather than pressure (although there are exceptions), a chief feature distinguishing dark energy from dark matter. (The constraint curves obtained in Ref. Schuecker:2002yj , e.g. Figure 3, are consistent with, but do not decisively pare down the parameter regions of our current results.) Furthermore, the analysis of mass power spectrum observations will have to take into account possible influences of a time-dependent $`w`$ and early quintessence, which we put off for later investigation. Overall, we have simulated in detail more than 400,000 individual cosmological models. The stored spectra and parameter-space likelihood functions will be used to evaluate additional constraints that can offer clues to the behavior of the dark energy. The set of thirteen target models, listed in Tables 1-4, which includes a revised concordance $`\mathrm{\Lambda }`$ model Wang:1999fa , should be of use by researchers to further test the quintessence hypothesis.

Acknowledgments This work was supported by NSF grant PHY-0099543 at Dartmouth. We thank Pier Stefano Corasaniti for useful conversations, and Dartmouth colleagues Barrett Rogers and Brian Chaboyer for use of computing resources.