text:
# A comparison of central temperatures of the intracluster gas determined from X-ray and SZ measurements

## 1 Introduction

While the X-ray spectroscopic measurement is a powerful and unique way nowadays to obtain the temperature of hot, diffuse gas contained within clusters of galaxies, it yields an emission-weighted temperature rather than the true temperature distribution. The lack of detailed information on the temperature profiles in clusters is probably the major source of uncertainties in the present determinations of the total masses and baryon fractions of clusters, which further hinders the dynamical properties of clusters from the cosmological applications such as the estimates of the cosmological parameters and the test of various models of structure formation in the universe. An independent, complementary method of measuring the distribution of intracluster gas is thus desirable. Silk & White (1978) were the first to suggest the utilization of the nonparametric reconstruction of the radial profiles of density \[$`n_e(r)`$\] and temperature \[$`T(r)`$\] of intracluster gas by combining the X-ray and SZ measurements. The basic idea is to inverse the observed X-ray and SZ temperature surface brightness profiles of clusters which have the functional forms of $`n_e^2(r)T^{1/2}(r)`$ and $`n_e(r)T(r)`$, respectively. And a simple combination of the two functions gives straightforwardly the gas density $`n_e(r)`$ and temperature $`T(r)`$. This method requires no assumption about the dynamical properties of clusters and the equation of state for intracluster gas, and can therefore be regarded as an ideal and ultimate tool to probe the gas distribution in clusters under spherical approximation.

Despite its elegant mathematical treatment further developed by Yoshikaya & Suto (1999) based on some theoretical models, the pioneering suggestion of Silk & White (1978) has not yet been put into practice to date. This has been primarily limited by the instrumental sensitivity and resolution of detecting the temperature variations of the cosmic background radiation (CBR) behind clusters. So far, the marginal detections of the radial SZ temperature distributions have been reported only for a few clusters (see Birkinshaw 1999), which can hardly be used for the purpose of reconstructing the gas temperature profiles because of the sparse data points and the large associated uncertainties. Nevertheless, we notice that the central measured or estimated SZ effects have been available for a great number of clusters with the past two decades’ efforts, and these central CBR temperature decrements in the Rayleigh-Jeans limit are usually used for setting constraints on the cosmological parameters ($`H_0`$, $`\mathrm{\Omega }_M`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$) in conjunction with the X-ray imaging and spectral measurements. Yet, another application of the central SZ data is instead to estimate the central temperature of intracluster gas when combined with the X-ray imaging observation, although this does not completely achieve the original goal of Silk & White (1978). It provides an indirect measurement of the X-ray temperature at cluster centers, which can be directly compared with the result given by the X-ray spectral analysis. Additionally, because both X-ray emission and SZ effect result from the gas distribution projected along the light of sight, a comparison of the central temperatures indicated by the two methods may also reveal valuable information about the radial temperature variations. Using the published X-ray and SZ data in literature, we present for the first time such a comparison in this paper and discuss its implications for dynamical properties of clusters. Throughout the paper we assume $`H_0=50`$ km s-1 Mpc-1 and $`\mathrm{\Omega }_0=1`$.

## 2 The model

The X-ray imaging observation provides a reliable measurement of the X-ray azimuthally-averaged surface brightness profile of a cluster which is usually described by the conventional $`\beta `$ model (Cavaliere & Fusco-Femiano 1976)

$$S_x(r)=S_0\left(1+\frac{r^2}{r_c^2}\right)^{3\beta +1/2}.$$
(1)
In the scenario of an optically-thin, thermal bremsstrahlung emission, the above form of $`S_x(r)`$ indicates (Cowie, Henriksen & Mushotzky 1987)

$$n_e(r)T^{1/4}(r)=n_{e0}T_0^{1/4}\left(1+\frac{r^2}{r_c^2}\right)^{3\beta /2},$$
(2)
where $`n_e`$ and $`T`$ are the electron number density and temperature, respectively. If we assume an equation of state, $`T(r)=T_0[n_e(r)/n_{e0}]^{\gamma 1}`$, we can write the electron number density as $`n_e(r)=n_{e0}(1+r^2/r_c^2)^\delta `$, in which $`\delta =6\beta /(3+\gamma )`$. The central electron number density $`n_{e0}`$, temperature $`T_0`$ and the X-ray surface brightness $`S_0`$ are connected by

$$n_{e0}^2=\left(\frac{3m_e\mathrm{}c^2}{2^4e^6}\right)\left(\frac{3m_ec^2}{2\pi kT_0}\right)^{1/2}\frac{4\pi ^{1/2}}{\mu _eg}\frac{\mathrm{\Gamma }(3\beta )}{\mathrm{\Gamma }(3\beta 1/2)}\frac{S_0(1+z)^4}{r_c},$$
(3)
where $`\mu _e=2/(1+X)`$, $`X=0.768`$ is the primordial hydrogen mass fraction, $`g1.2`$ is the average Gaunt factor, $`\mathrm{\Gamma }`$ is the gamma function, and $`z`$ is the cluster redshift.

On the other hand, the CBR temperature decrement/increment at the cluster center predicted by the SZ effect is

$`{\displaystyle \frac{\mathrm{\Delta }T_{sz}}{T_{\mathrm{CBR}}}}=4\xi (x){\displaystyle _0^+\mathrm{}}{\displaystyle \frac{kT(r)}{m_ec^2}}\sigma _Tn_e(r)𝑑r,`$ (4)
$`\xi (x)={\displaystyle \frac{x^2e^x}{2(e^x1)^2}}\left(4x\mathrm{coth}{\displaystyle \frac{x}{2}}\right),`$ (5)
where $`T_{\mathrm{CBR}}`$ is the temperature of the present CBR, and $`x=h\nu /kT_{\mathrm{CBR}}`$ is the dimensionless frequency. A straightforward computation using eqs.(2)-(5) as well as the polytropic equation of state gives

$`\left({\displaystyle \frac{kT_0}{m_ec^2}}\right)^{3/2}=`$ $`\left({\displaystyle \frac{\mathrm{\Delta }T_{sz}}{T_{\mathrm{CBR}}}}\right)^2\left({\displaystyle \frac{\mathrm{\Gamma }(3\beta ^{}/2)}{\mathrm{\Gamma }(3\beta ^{}/21/2)}}\right)^2{\displaystyle \frac{\mathrm{\Gamma }(3\beta 1/2)}{\mathrm{\Gamma }(3\beta )}}`$ (6)
$`{\displaystyle \frac{\sqrt{2}\alpha \mu _egm_ec^3}{8\pi ^2\sqrt{3}(3\beta 3/2)\xi ^2\sigma _T(1+z)^4S_0r_c}}`$
in which $`\beta ^{}=3\beta \gamma /(3+\gamma )`$, and $`\alpha `$ is the fine structure constant. The above expression is often used in the determination of the Hubble constant for an isothermal gas distribution if we write $`r_c=d_A(z)\theta _c`$ where $`d_A`$ is the angular diameter distance to the cluster at $`z`$. On the other hand, given the Hubble constant and the polytropic index $`\gamma `$, all the parameters at the right-hand of eq.(6) are measurable from the X-ray imaging and SZ measurements, which will instead allow us to derive the central temperature of the X-ray emitting gas.

## 3 Application to X-ray clusters

Our cluster sample consists primarily of 21 nearby clusters studied extensively by Mason & Myers (2000) plus 10 high-redshift clusters with recent, reliable SZ data in the literature. Inclusion of a cluster is also subject to whether or not the good spectral temperature data are available in order to facilitate our comparison between different methods. We take the best-fit values of $`\beta `$, $`r_c`$ and $`S_0`$ (or $`L_x`$) as well as the emission-weighted temperature data directly from the literature. For the latter, we adopt the temperature given by a single-phase (s.p.) model and the value by excluding the cooling flows (c.l.) (White 2000), respectively. Since the X-ray spectrum is always observed over a finite bandpass, we assume an optically-thin, thermal bremsstrahlung emission model with the metal abundances of $`30\%`$ solar to convert the observed $`S_0`$ (or $`L_x`$) over the corresponding energy band into the bolometric value. We use the central SZ decrement/increment of $`\mathrm{\Delta }T_{sz}`$ after the correction of the finite beam in each observation. All the input data of the X-ray and SZ measurements are summarized in Table 1.

For an isothermal gas distribution, the jointly determined central temperature, $`T_{est}`$, from the X-ray imaging and SZ measurements (eq.) can be directly compared with the X-ray spectral result, $`T_{spec}`$. We display in Fig.1 our derived $`T_{est}`$ versus the observed $`T_{spec}`$ (s.p. model) for the 31 clusters in Table 1. A glimpse of Fig.1 reveals that the two sets of data are roughly consistent with each other. The best-fit $`T_{spec}`$-$`T_{est}`$ relation reads

$`T_{est}=10^{0.22\pm 0.18}T_{spec}^{1.22\pm 0.20},`$ $`\mathrm{s}.\mathrm{p}.\mathrm{model};`$ (7)
$`T_{set}=10^{0.36\pm 0.26}T_{spec}^{1.29\pm 0.29},`$ $`\mathrm{c}.\mathrm{l}.\mathrm{corrected},`$ (8)
in which the error bars are the $`90\%`$ confidence limits which have taken the measurement uncertainties into account. We have also tested the dependence of the best-fit $`T_{spec}`$-$`T_{est}`$ relation on the Hubble constant by adopting a larger value of $`H_0=75`$ km s-1 Mpc-1. This only leads to a minor modification to the above results. For instance, the best-fit relation for the s.p. model now reads $`T_{est}=10^{0.34\pm 0.11}T_{spec}^{1.22\pm 0.12}`$.

If the intracluster gas deviates from isothermality, our comparison of $`T_{est}`$ and $`T_{spec}`$ becomes to be a little complex because we need to recover the central temperature from the emission-weighted spectral fit. On the other hand, we can place a useful constraint on the polytropic index $`\gamma `$ if we demand that our derived central temperature from eq.(6) identify the spectral value after a proper correction to the emission-weighted temperature $`\overline{T}`$ predicted by

$$\frac{dE}{d\nu }\overline{T}^{1/2}g_{ff}(\overline{T},\nu )\mathrm{exp}(h\nu /k\overline{T})n_e^2𝑑V,$$
(9)
where $`g_{ff}`$ is the Gaunt factor of the free-free emission, $`\nu `$ is the X-ray observing frequency, and the integration is performed over the whole cluster volume if we neglect the finite detection aperture. This compares with the theoretically expected spectrum for the hot intracluster gas

$$\frac{dE}{d\nu }n_e^2T^{1/2}g_{ff}(T,\nu )\mathrm{exp}(h\nu /kT)𝑑V.$$
(10)
For each cluster with the reported $`\overline{T}=T_{spec}`$, we use the Monte Carlo simulations over the typical energy band $`0.5`$ \- $`10`$ keV to obtain the central temperature $`T_0`$ and the corresponding $`\gamma `$ with the restriction of eq.(6). The resultant $`\gamma `$ value for each cluster is listed in Table 1, and a combined analysis of the results among the 31 clusters yields

$`\gamma =0.91\pm 0.14,`$ $`\mathrm{s}.\mathrm{p}.\mathrm{model};`$ (11)
$`\gamma =0.93\pm 0.14,`$ $`\mathrm{c}.\mathrm{l}.\mathrm{corrected},`$ (12)
which indicates that clusters are essentially consistent with isothermality at $`90\%`$ significance level.

## 4 Discussion and conclusions

Our joint analysis of the indirect and direct temperature measurements from the X-ray and SZ observations provides a useful clue to resolving the temperature profile discrepancy raised in recent years particularly by Markevitch et al. (1998), Irwin, Bregman & Evrard (1999), White (2000) and Irwin & Bregman (2000). Markevitch et al. (1998) claimed a significant temperature decline with increasing radius represented by a polytropic index of 1.2-1.3, while the other investigators have essentially detected an isothermal temperature profile. Our present study gives an independent yet convincing evidence for the assumption of isothermality. The advantage of our algorithm over the previous analysis is that we have obtained the central temperatures of clusters without utilizing the X-ray spectral data, while a combination of both has allowed us to set a valuable constraint on the radial temperature profiles of the hot X-ray emitting gas.

A joint determination of the central temperature of intracluster gas from the X-ray imaging and SZ measurements can be used at present as an independent and complementary method to examine the reliability and accuracy of the temperature measurement provided by the conventional X-ray spectral analysis, although our ultimate goal is to obtain the detailed radial profiles of electron number density and temperature when the SZ surface brightness can be precisely mapped (Silk & White, 1978; Yoshikawa & Suto 1999). Yet, we should point out that it may not be easy to identify the temperature structure of the hot gas from such a combined analysis, in the sense that the dominant component in both $`S_x`$ and $`\mathrm{\Delta }T_{sz}`$ is characterized by the sharp decrease of the electron number density with outward radius, while it is unlikely that the temperature profiles can demonstrate a similar variation. An accurate measurement of $`S_x`$ and $`\mathrm{\Delta }T_{sz}`$ at outermost radii will thus be needed if one wants to apply the present technique to the reconstruction of temperature profiles of intracluster gas. Alternatively, the present technique requires that the X-ray surface brightness is well represented by the $`\beta `$ model. Presence of cooling flows in the central regions of many X-ray clusters may yield an underestimate of the $`\beta `$ parameter and the core radius $`r_c`$ and an overestimate of the central surface brightness $`S_0`$, which will in turn affect our estimate of $`T_{est}`$ through eq.(6). Indeed, an examination of Fig.1 reveals that all the massive cooling flow clusters with $`\dot{M}>600`$ $`M_{}`$ yr-1 (A478, A1835, A2204 and Zw3146) demonstrate rather large scatters around the $`T_{est}=T_{spec}`$ relation. Therefore, the present method is most applicable to morphologically simple clusters without cooling flows.

This work was supported by the National Science Foundation of China, under Grant 19725311.


format: text

text:
# Nonlinear Stiffness, Lyapunov Exponents, and Attractor Dimension

(1999)

## Abstract

I propose that stiffness may be defined and quantified for nonlinear systems using Lyapunov exponents, and demonstrate the relationship that exists between stiffness and the fractal dimension of a strange attractor: that stiff chaos is thin chaos.

journal: Physics Letters Avolume: 264

,

What constitutes a stiff dynamical system? Stiffness is closely related to numerical methods : the signature of stiffness in a problem is that, upon integration with a general numerical scheme — a method not specially designed for stiff problems — the routine takes extremely small integration steps , which makes the process computationally expensive. One view is that stiffness is inextricably linked with the numerical integration scheme used, so that there would be no such thing as an intrinsically stiff dynamical system, and the best we could hope for is an operational definition such as that above . Moreover, it has been proposed that chaotic problems cannot be stiff . I argue below that this is not the case, and provide a definition and a quantitative measure of stiffness for nonlinear dynamical systems. I demonstrate how stiffness affects the geometry of the strange attractor of a chaotic system: that stiff chaos is thinner — has smaller fractional part of the fractal dimension — than nonstiff chaos.

When integrating a stiff problem with a variable-step explicit numerical integration scheme, the initial step length chosen causes the method to be at or near numerical instability, which leads to a large local truncation error estimate. This causes the numerical routine to reduce the step length substantially, until the principal local truncation error is brought back within its prescribed bound. The routine then integrates the problem successfully, but uses a far greater number of steps than seems reasonable, given the smoothness of the solution. Because of this, computation time and round-off error are a problem when using conventional numerical integration techniques on stiff problems, and special methods have been developed for them.

Traditionally in numerical analysis, a linear stiff system of size $`n`$ is defined by

$$\mathrm{Re}(\lambda _i)<0,1in,$$
(1)
where $`\lambda _i`$ are the eigenvalues of the Jacobian of the system, with

$$\underset{1in}{\mathrm{max}}|\mathrm{Re}(\lambda _i)|\underset{1in}{\mathrm{min}}|\mathrm{Re}(\lambda _i)|.$$
(2)
The stiffness ratio $`R`$ provides a quantitative measure of stiffness:

$$R=\frac{\underset{1in}{\mathrm{max}}|\mathrm{Re}(\lambda _i)|}{\underset{1in}{\mathrm{min}}|\mathrm{Re}(\lambda _i)|}.$$
(3)
By this definition, a stiff problem has a stable fixed point with eigenvalues of greatly differing magnitudes; large negative eigenvalues correspond to fast-decaying transients $`e^{\lambda t}`$ in the solution.

As an example of a linear stiff problem, consider the equation

$$y^{\prime \prime }+1001y^{}+1000y=0.$$
(4)
We can write this as the vector equation $`𝐲^{}=A𝐲`$ where

$$A=\left(\begin{array}{cc}0& 1\\ 1000& 1001\end{array}\right),$$
(5)
and the eigenvalues are $`\lambda _1=1`$ and $`\lambda _2=1000`$. This has solution

$$y=Ae^t+Be^{1000t},$$
(6)
so when integrating the problem numerically with a general variable-step method we would expect to be able to use a large integration step length after the $`e^{1000t}`$ transient term decays, but in fact the presence of the large negative eigenvalue $`\lambda _2`$ prevents this. With appropriate initial conditions, one can even set $`B=0`$ and remove the $`e^{1000t}`$ term from the solution entirely; one nevertheless has to use a very small step length throughout the calculation, as step length is still dictated by the size of $`\lambda _2`$.

The definition of linear stiffness of Eqs. (2) and (3) is not relevant for nonlinear systems. The stiffness ratio Eq. (3) is often not a good measure of stiffness even for linear systems, since if the minimum eigenvalue is zero, the problem has infinite stiffness ratio, but may not really be stiff at all if the other eigenvalues are of moderate size. The inadequacy of the stiffness ratio is clearly recognized by numerical analysts, who have moved away from trying to pin down the definition of stiffness and have generally adopted the pragmatic approach I alluded to earlier: “Stiff equations are problems for which explicit methods don’t work” .

Let us look at a nonlinear stiff problem. In Fig. 1, I show the results of integrating the van der Pol equation

$$y^{\prime \prime }\mu (1y^2)y^{}+y=0$$
(7)
using a standard variable-step fourth-order Runge–Kutta code , for the cases $`\mu =1`$ and $`\mu =100`$. There is the same bound on the principal local truncation error estimate in both cases. We can see from the far greater number of steps needed for $`\mu =100`$, that at large $`\mu `$ the van der Pol equation becomes very stiff. The steps are so small at $`\mu =100`$ that the individual crosses representing each step merge to form a continuous broad line on the graph. At these large $`\mu `$ values, the equation describes a relaxation oscillator . These have fast and slow states in their cycle which characterizes the jerky motion displayed in Fig. 1(b). And if we introduce forcing into the van der Pol equation

$$y^{\prime \prime }\mu (1y^2)y^{}+y=A\mathrm{cos}\omega t$$
(8)
we obtain a chaotic system  which, like its unforced counterpart, is stiff, as demonstrated in Fig. 1(c)), where the steps are so small that they have merged together in the image.

Throughout the foregoing, I have illustrated stiff systems with reference to their effect on the numerical method used for integrating them. Is it possible then to define nonlinear stiffness independently of the numerical method used? Lambert  puts forward the definition: A system is said to be stiff in a given interval $`x`$ if in that interval the neighbouring solution curves approach the solution curve at a rate which is very large in comparison with the rate at which the solution varies in that interval. This definition does not make reference to the numerical method used and is instead concerned with the curvature and local Lyapunov exponents of the solution. The curvature of the solution curve

$$\kappa =\frac{y^{\prime \prime }}{(1+y^2)^{3/2}}$$
(9)
at a point quantifies the wiggliness of the trajectory at that point, while local Lyapunov exponents are defined as

$$\gamma _i(\tau ,t)=\underset{\sigma _i(\tau )0}{lim}\frac{1}{\tau }\mathrm{ln}\frac{\sigma _i(t+\tau )}{\sigma _i(t)},$$
(10)
where $`\sigma _i(t)`$ are the principal axes of an ellipsoidal ball evolving in time in phase space. The local Lyapunov exponent $`\gamma _i(\tau ,t)`$ is a generalization of the Lyapunov exponent $`\lambda _i`$: the Lyapunov exponent is the limit of $`\tau `$ going to infinity in the local Lyapunov exponent

$$\lambda _i=\underset{\tau \mathrm{}}{lim}\gamma _i(\tau ,t).$$
(11)
Lyapunov exponents show the rate of convergence or divergence of neighbouring trajectories, and an $`n`$-dimensional system has $`n`$ Lyapunov exponents corresponding to $`n`$ expanding or contracting directions in phase space. Since Lyapunov exponents are defined in the infinite time limit, they cannot reflect differing rates of convergence or divergence in different parts of a trajectory. Whereas the Lyapunov exponent is the same for almost all starting points on a trajectory, the local Lyapunov exponent can vary depending on the starting point and the length of trajectory examined. Fast convergence to a neighbouring trajectory, which implies having large negative local Lyapunov exponents, indicates stiffness. A system that is equally stiff at all points along a trajectory has a constant convergence rate, and in this case, the local Lyapunov exponent will be the same as the Lyapunov exponent. Often, though, a system can show intervals of stiff and nonstiff behaviour. Relaxation oscillators like the van der Pol oscillator are examples of this, having a stiff slow manifold, and a nonstiff fast manifold. With each oscillation we have two intervals of stiff behaviour interspersed with two intervals of nonstiff behaviour.

As it stands, Lambert’s definition of stiffness will not do for us, since it assumes that the system is not chaotic: the principal local Lyapunov exponent is large and negative to obtain fast convergence of neighbouring trajectories. Instead, we need to adapt it to allow for chaotic behaviour along with stiffness, by looking at the largest negative local Lyapunov exponent, rather than the principal local Lyapunov exponent. Our definition of nonlinear stiffness is then: A system is stiff in a given interval if in that interval the most negative local Lyapunov exponent is large, while the curvature of the trajectory is small. A quantitative measure of nonlinear stiffness at any point can then be obtained from the ratio of the most negative local Lyapunov exponent and the curvature of the trajectory:

$$R_{nl}=\frac{\left|\underset{1in}{\mathrm{min}}\gamma _i(\tau ,t)\right|}{\kappa (t)}.$$
(12)
If desired this could be averaged over the trajectory to give a measure of mean stiffness.

We now have a definition of nonlinear stiffness; what does it imply for chaotic systems? If a Lyapunov exponent of a system is large and negative, then the local Lyapunov exponent must be large and negative at least over some of the trajectory, so, given suitable bounds on the curvature of the trajectory, a large negative Lyapunov exponent is a sufficient condition for stiffness. On the other hand, at least one positive Lyapunov exponent is necessary for chaos. Thus for stiff chaos we should have a large spread of Lyapunov exponents, with at least one positive, and one large and negative.

Lyapunov exponents are certainly related to the fractal dimension of an attractor. The Kaplan–Yorke conjecture  holds that the fractal dimension is

$$D_{KY}=\frac{1}{|\lambda _{j+1}|}\underset{k=1}{\overset{j}{}}\lambda _k+j,$$
(13)
where $`j`$ is the expansion dimension of the system: the largest integer such that

$$\underset{k=1}{\overset{j}{}}\lambda _k0.$$
(14)
The Kaplan–Yorke estimate is often good — i.e., very close to the measured (capacity) dimension. In general the fractal dimension $`D`$ lies between the expansion dimension and the Kaplan–Yorke estimate :

$$jDD_{KY}.$$
(15)
For example, in a three-dimensional chaotic system, $`\lambda _1>0`$, $`\lambda _2=0`$, and $`\lambda _3<0`$, so

$$D_{KY}=2+\frac{\lambda _1}{|\lambda _3|}.$$
(16)
From this, if $`|\lambda _3|`$ gets larger, with increasing stiffness, then the fractal dimension of the strange attractor moves closer to two.

The size of the fractional part of the fractal dimension of a strange attractor may be termed the thickness or thinness of the chaos. The folded structure of the foliations of the strange attractor — the Abraham and Shaw bagel  — is tightly wound in the case of thin chaos with small fractional part. On the other hand, if the fractal folding of the surface of the strange attractor is macroscopic, the chaos is thick chaos having a large fractional part to its fractal dimension.

We should expect then that stiff chaos will be thin, and this indeed proves to be the case. In Fig. 1 I illustrated the stiffness of the forced van der Pol equation. This equation is also well known for the difficulty of capturing pictorially its chaotic nature at large $`\mu `$ . On the other hand, with the Shaw variant of the forcing , which may be written as a three-dimensional autonomous system

$`\dot{x}`$ $`=`$ $`y+\mu (1y^2)x`$
$`\dot{y}`$ $`=`$ $`x+A\mathrm{cos}\omega z`$
$`\dot{z}`$ $`=`$ $`1`$ (17)
whose divergence — the sum of the Lyapunov exponents — is $`\mu (1y^2)`$, a strange attractor is easily found at $`A=0.932`$, $`\mu =1.18`$, and $`\omega =1.86`$. Now, with the benefit of the above analysis, we see how fractal dimension, dissipation, and stiffness are all linked together: whereas the normal forced van der Pol equation studied at large $`\mu `$ (high dissipation) is stiff, making the chaos thin, the attractor in the Shaw variant at small $`\mu `$ (low dissipation) displays thick chaos, and, following these arguments, will not be stiff; this is indeed the case (Fig. 2).

I should like to thank the referees for their helpful comments. I acknowledge the financial support of the Spanish Consejo Superior de Investigaciones Científicas.


text:
# References

A Quantum Fluctuation Theorem

Jorge Kurchan

P.M.M.H. Ecole Supérieure de Physique et Chimie
Industrielles, 10 rue Vauquelin 75231 Paris, France

## Abstract

We consider a quantum system strongly driven by forces that are periodic in time. The theorem concerns the probability $`P(e)`$ of observing a given energy change $`e`$ after a number of cycles. If the system is thermostated by a (quantum) thermal bath, $`e`$ is the total amount of energy transferred to the bath, while for an isolated system $`e`$ is the increase in energy of the system itself. Then, we show that $`P(e)/P(e)=e^{\beta e}`$, a parameter-free, model-independent relation.

In the past few years there has been a renewed interest in the study of quantum systems out of equilibrium, to a large extent stimulated by the design of new experimental settings and by the construction of new devices. If a system is well out of equilibrium, as for example when it is strongly driven by periodic forces, then linear response theory (understood as linear perturbations around the Gibbs measure) is insufficient.

Even in the context of classical mechanics not many generic results are available beyond linear response. An interesting new development consists of a number of relations for strongly out of equilibrium systems, mainly regarding the distribution of work and entropy production. The first of such fluctuation theorems was discovered by Evans et.al , who understood that the basic ingredient was time-reversal symmetry.

Two important further steps, made by Gallavotti and Cohen widened the scope and interest of the subject. On the one hand, it was realized  that the fluctuation theorems are indeed the far from equilibrium generalisations of the well-known equilibrium theorems (fluctuation-dissipation and Onsager reciprocity). Most intriguingly, a byproduct of their proof  was that, just as the validity of the fluctuation-dissipation relation is a strong indication of equilibration, the fact that a fluctuation formula holds in a driven stationary system strongly hints that the system can be considered ‘as ergodic as possible’ — with all the implications this entails.

These results concern deterministic systems. If instead a finite system is in contact with a stochastic thermal bath, then the ‘ergodicity’ questions become trivial, and all the fluctuation formulae are extremely simple to prove .

Particularly relevant for the present work are the simple and remarkable Jarzynski and other ‘work relations’ valid well out of equilibrium , whose close relation to the fluctuation theorems was clarified by Crooks .

Except for Ref., the developments described so far are restricted to classical mechanics. The purpose of this paper is to prove two versions of a fluctuation theorem for quantum systems under strong periodic drive, either isolated or in contact with a thermal bath. Very recently, these questions have become relevant in the context of detection of quantum ‘shot’ noise generated by currents flowing through devices, given the possibility it offers to observe different quasiparticle charges of the carriers .

We shall consider an evolution generated by a time-dependent Hamiltonian $`H(q,p,t)`$ with real matrix elements, so that there is time-reversal symmetry in each infinitesimal time step. We shall assume the explicit time-dependence is periodic (with period $`t_o`$) and that the cycles are symmetric :

$$H(q,p,t)=H(q,p,t)=H(q,p,t)$$
(1)
Splitting the evolution in infinitesimal steps $`t_1,\mathrm{}t_n`$:

$$|\varphi (t_f)=U|\varphi (t_i)=U_nU_{n1}\mathrm{}U_2U_1|\varphi (t_i)$$
(2)
where $`n`$ is large, $`U_re^{\frac{i}{\mathrm{}}\tau H_r}`$ and $`\tau =(t_ft_i)/n0`$, the reality condition implies that: $`U_r^{}=U_r^{}`$ and, together with the symmetry of cycles Eq. (1) this implies that the evolution over a cycle satisfies time-reversibility: Consider the case in which the time $`t_ft_i`$ consists of an integer number of cycles, and hence $`U_r=U_{nr+1}`$ (cfr. (1)). Then:

$`U^{}`$ $`=`$ $`U_1^{}U_2^{}\mathrm{}U_{n1}^{}U_n^{}=\left[U_1U_2\mathrm{}U_{n1}U_n\right]^{}`$ (3)
$`=`$ $`\left[U_nU_{n1}\mathrm{}U_2U_1\right]^{}=U^{}`$
a formula only valid for $`(t_ft_i)/t_o=integer`$, which we shall assume throughout this paper.

We denote the eigenvectors of the initial time Hamiltonian $`|\psi _\alpha `$:

$$H_1|\psi _\alpha =\epsilon _\alpha |\psi _\alpha $$
(4)
and the corresponding partition function:

$$Z=\text{Tr}e^{\beta H_1}$$
(5)
Consider the following protocol:

* With (canonical) probability $`p_\alpha =e^{\beta \epsilon _\alpha }/Z`$ we choose a wavefunction $`|\psi _\alpha `$.
* We let it evolve through an integer number of cycles:

$$|\varphi =U|\psi _\alpha =\underset{\gamma }{}|\psi _\gamma \psi _\gamma |\varphi $$
(6)
* We measure the final $`H_1`$ and record the energy difference $`e`$ between initial and final times.

Let us calculate the probability distribution of $`e`$. For a given value of $`\alpha `$, this distribution reads:

$`P_\alpha (e)`$ $`=`$ $`{\displaystyle \underset{\gamma }{}}\delta [e(\epsilon _\gamma \epsilon _\alpha )]\left|\psi _\gamma |\varphi \right|^2`$ (7)
$`=`$ $`{\displaystyle \underset{\gamma }{}}\delta [e(\epsilon _\gamma \epsilon _\alpha )]\left|\psi _\gamma |U|\psi _\alpha \right|^2`$
The average distribution $`P(e)`$ over initial conditions is:

$$P(e)=\frac{1}{Z}\underset{\gamma ,\alpha }{}\delta [e(\epsilon _\gamma \epsilon _\alpha )]\left|\psi _\gamma |U|\psi _\alpha \right|^2e^{\beta \epsilon _\alpha }$$
(8)
Writing the delta function in integral form, we get:

$`P(e)`$ $`=`$ $`{\displaystyle \frac{1}{Z}}{\displaystyle _i\mathrm{}^i\mathrm{}}𝑑\lambda e^{\lambda [e(\epsilon _\gamma \epsilon _\alpha )]}{\displaystyle \underset{\gamma ,\alpha }{}}\left|\psi _\gamma |U|\psi _\alpha \right|^2e^{\beta \epsilon _\alpha }`$ (9)
$`=`$ $`{\displaystyle \frac{1}{Z}}{\displaystyle _i\mathrm{}^i\mathrm{}}𝑑\lambda e^{\lambda e}Q(\lambda )`$
where:

$`Q(\lambda )`$ $`=`$ $`{\displaystyle \underset{\gamma ,\alpha }{}}e^{\lambda (\epsilon _\gamma \epsilon _\alpha )}\left|\psi _\gamma |U|\psi _\alpha \right|^2e^{\beta \epsilon _\alpha }`$ (10)
$`=`$ $`{\displaystyle \underset{\gamma ,\alpha }{}}e^{\lambda \epsilon _\gamma }\psi _\alpha |U^{}|\psi _\gamma \psi _\gamma |U|\psi _\alpha e^{(\lambda +\beta )(\epsilon _\alpha )}`$
$`=`$ $`\text{Tr}\left[U^{}e^{\lambda H_1}Ue^{(\lambda +\beta )H_1}\right]`$
The fluctuation theorem can be proved very simply. Firstly, let us prove the following KMS-like  relation. Let $`D`$ be any real symmetric operator. Then:

$`Q_D(\lambda )`$ $``$ $`\text{Tr}\left[U^{}e^{\lambda D}Ue^{(\lambda +\beta )D}\right]=\text{Tr}\left[U^{}e^{(\lambda +\beta )D}Ue^{\lambda D}\right]`$ (11)
$`=`$ $`Q_D(\lambda ^{}\beta )^{}=Q_D(\lambda \beta )`$
where we have used transposition, cyclic permutation and Eq. (3). Putting $`D=H_1`$, this is a form of the fluctuation theorem for $`Q(\lambda )=Q_D(\lambda )`$.

In order to see what the implications of (11) are for $`P(e)`$ we shall need to show that $`Q(\lambda )`$ is analytic on the stripe $`\beta \text{Re}(\lambda )0`$. To do this, we use the fact that:

$$\left|\text{Tr}AB\right|^22\left(\text{Tr}AA^{}\right)\left(\text{Tr}BB^{}\right)$$
(12)
Putting $`A=U^{}e^{\lambda H_1}`$ and $`B=Ue^{(\lambda +\beta )H_1}`$:

$$\left|Q(\lambda )\right|^22\left(\text{Tr}e^{2\left\{\text{Re}(\lambda )+\beta \right\}H_1}\right)\left(\text{Tr}e^{2\text{Re}(\lambda )H_1}\right)$$
(13)
Because the partition function (5) converges for positive temperatures, neither factor diverges if $`\beta <\text{Re}(\lambda )<0`$

Inserting (11) with $`D=H_1`$ in (9), the analyticity result allows us to shift the integration from $`\text{Re}(\lambda )=0`$ to $`\text{Re}(\lambda )=\beta `$, and we get:

$$P(e)=P(e)e^{\beta e}$$
(14)
This is the fluctuation theorem for the probability of energy changes in an isolated system. Note that we could have obtained this result directly, without writing the KMS equation, by using the (time-reversal) symmetry of $`U`$ (equation (3)) in (8).

Taking the expectation value of $`e^{\beta e}`$ over an integer number of cycles, one obtains a quantum version of a Jarzynski work formula :

$$\overline{e^{\beta e}}=𝑑eP(e)e^{\beta e}=𝑑eP(e)=1$$
(15)
where the overline denotes average over quantum amplitudes and initial conditions. This is a rather surprising model-independent result, which we can rewrite as:

$$0=\frac{1}{\beta }\mathrm{ln}\overline{e^{\beta e}}\frac{1}{\beta }\overline{\mathrm{ln}e^{\beta e}}=\overline{e}$$
(16)
We obtain the second principle, arising as the familiar inequality

$$annealedaveragequenchedaverage$$
(17)
with the initial conditions playing the role of disorder.

Generalisation and Thermostated systems

In the most physical setting, we have a system in contact with a bath. The bath can be modeled for example with an infinite set of harmonic oscillators coupled to each variable in the system. Consider for simplicity the case of one system variable:

$$H(x,p_x,y_1,p_1,\mathrm{},y_M,p_M)=H_{system}+H_{int}+H_{bath}$$
(18)
with:

$$H_{bath}=\underset{i}{}\left[\frac{p_i^2}{2m}+m\omega _i^2\frac{y_i^2}{2}\right]$$
(19)
$$H_{int}=\underset{i}{}C_iy_ix$$
(20)
and, say:

$$H_{system}=\frac{p_x^2}{2m_x}+V(x,t)$$
(21)
(The interested reader will find an extensive literature on this implementation of heat baths in  and references therein.)

We can define the following bases of eigenvalues:

$`H_{bath}|\chi _\alpha `$ $`=`$ $`\epsilon _\alpha ^y|\chi _\alpha `$
$`O_{system}|\psi _\alpha ^{}`$ $`=`$ $`o_\alpha ^{}|\psi _\alpha ^{}`$ (22)
$`O`$ is any real symmetric operator corresponding to an observable depending exclusively on the system variables, with a spectrum bounded from below (e.g. $`H_{system}(t_i)`$). The wavefunctions $`|\chi _\alpha `$ and $`|\psi _\alpha `$ depend only on bath and system variables, respectively.

We consider an initial condition with no correlations between bath and system constructed as follows:

* We choose a bath wavefunction $`|\chi _\alpha `$ with canonical probability $`e^{\beta \epsilon _\alpha ^y}`$.
* We choose a system wavefunction $`|\psi _\alpha ^{}`$ with probability $`p_\alpha ^{}e^{\beta o_\alpha ^{}}`$. This distribution is quite general, given the freedom of choice of $`O`$. It could be for example a canonical distribution at higher temperature ($`O=\frac{\beta ^{}}{\beta }H_1`$), as resulting from a temperature quench.
* We start with the initial state $`|\chi _\alpha |\psi _\alpha ^{}`$ and let it evolve.
* We measure

$$DH_{bath}+O$$
(23)
at the beginning and at the end, and record the difference. Note that at this stage the operator $`O`$ which we are measuring is forced by the initial condition we are choosing.

It is now easy to prove that the fluctuation formula (14) holds, with $`e`$ measuring the difference in value of bath energy plus observed value of $`O`$. To do this one uses (11) with $`D`$ as in (23). In order to guarantee the analyticity of $`Q_D`$, one must assure that (see argument leading to (13))

$`\text{Tr}\left[e^{2\left\{\text{Re}(\lambda )+\beta \right\}D}\right]`$ $`=`$ $`\text{Tr}\left[e^{2\left\{\text{Re}(\lambda )+\beta \right\}H_{bath}}\right]`$ (24)
$`\times `$ $`\text{Tr}\left[e^{2\left\{\text{Re}(\lambda )+\beta \right\}O}\right]<\mathrm{}`$
and

$$\text{Tr}\left[e^{2\text{Re}(\lambda )D}\right]=\text{Tr}\left[e^{2\text{Re}(\lambda )H_{bath}}\right]\text{Tr}\left[e^{2\text{Re}(\lambda )O}\right]<\mathrm{}$$
(25)
for $`\beta <\text{Re}(\lambda )<0`$. The traces over the bath variables are bounded, since they correspond to the bath partition function. For the initial probabilities of the system, these conditions require that:

$$\text{Tr}\left[e^{2\beta \mu O}\right]\underset{\alpha }{}p_\alpha ^{2\mu }$$
(26)
is finite for $`0<\mu <1`$, a condition we assume.

We are now in a position of discussing a general system in contact with a an infinite heat bath after many cycles have elapsed ($`t_ft_i\mathrm{}`$). If we have that, under these circumstances:

i) The expected value of $`O`$ — a property of the system — stays finite as $`t_f`$ grows.

ii) The bath receives by virtue of the time-dependent forces an energy proportional to the number of cycles : $`e=(t_ft_i)e_o+e_1`$, with $`e_1`$ finite and and $`e_o`$ independent of the initial configuration.

Then, the derivation above of the fluctuation theorem carries over for long times to any initial distribution and to leading order in $`t_ft_i`$, equation (14) will be a statement about the probability distribution of the energy $`e_o`$ the bath received, i.e. its entropy increase . Measuring the fluctuation formula in a concrete situation becomes then a test for a property of the ‘stationary’ (i.e. periodical) asymptotic quantum state.

Acknowledgements

I wish to thank L.F. Cugliandolo, D. Grempel and P. Leboeuf for useful suggestions. I am indebted to Y. Imry for calling my attention to the connections with works on quantum devices.


text:
# Cosmic Acceleration: A Natural Remedy For Horizon and Flatness Problems

## I Introduction

Only a few years ago, it was believed that the universe is undergoing a declerated expansion (DE). Then it was felt that there is no compelling reason to retain the “cosmological constant” term, $`\mathrm{\Lambda }`$, in the Einstein equations describing the evolution of the universe. In such a case, it followed that, the behaviour of the universe filled with matter and radiation, could be broadly described as $`S(t)=\alpha t^n`$, where $`n<1`$ and $`\alpha `$ is a constant of proportionality. For instance, for the early radiation dominated era, one has $`St^{1/2}`$ and for the matter dominated universe, $`St^{2/3}`$ (for a closed universe with curvature parameter $`k=0`$).

However, over the past few years, studies of luminosity-distance studies involving type supernova I, have suggested that the universe is undergoing accelerated expansion (AE). Such AE is possible either due to presence of a $`\mathrm{\Lambda }`$ term or because of the presence of a quintessence. The simplest way to broadly reprsent the AE would be to write

$$S(t)=\alpha t^n;n>1$$
(1)
In the past era of supposed decelerated expansion, one would encounter the socalled horizon and flatness puzzle. In the former case, it is difficult to find an answer asto why the universe is so isotropic and homogeneous if the universe comprised causally disconnected tiny patches in the very early era. In the latter case, one has to explain why the quantity $`\mathrm{\Omega }1𝒪1`$ in the present epoch unless $`\mathrm{\Omega }`$ was extremely fine tuned to unity in the early past. Here $`\mathrm{\Omega }`$ is the ratio of the actual density of the universe ($`u`$ to the critical density needed to close the universe$`u_c=8\pi G/3H^2`$, where $`H`$ is the Hubble constant. To remedy such problems, in the framework of the supposed (incorrect) decelerated expansion, one needed to introduce, in a fairly ad hoc manner, a brief spell of sudden exponential expansion by a factor of $`10^{30}`$ starting from $`t10^{35}`$ s. It is this rapid expansion which is supposed to to bring distant regions within the sphere of causality and simultaneously flattening the universe. Following this brief infalationary period, the universe is supposed to revert back to the original decelerated expansion mode in a manner which is harly well understood.

However, we shall show below that once we accept the fact the universe is undergoing AE rather than DE and assume that it was so in the past too, there is no need to posulate an ad hoc inflationary phase alteast for the purpose of solving the horizon or flatness problems.

## II Removal of Horizon Problem

The range of the causal horizon is defined as

$$d_H=S(t)_0^t\frac{dt^{}}{S(t^{})}$$
(2)
Using Eq. (1) we see that

$$d_H=\frac{S(t)}{k}_0^tt^n𝑑t^{}=\frac{S(t)}{k}\frac{t^{(1n)}}{(1n)}_0^t$$
(3)
For $`n>1`$, the above integral blows up and there is no horzon problem.

## III Flatness Problem

In the presence of the cosmological term, the relationship between present density $`u`$, Hubble constant $`H`$, and scale size $`S`$ is

$$H^2+\frac{K}{S^2}=\frac{8\pi Gu}{3}+\frac{\mathrm{\Lambda }c^2}{3}$$
(4)
Dividing both sides by $`H^2`$, we have

$$1+\frac{K}{H^2S^2}=\mathrm{\Omega }+\frac{\mathrm{\Lambda }c^2}{3H^2}$$
(5)
Or,

$$\mathrm{\Omega }1=\frac{k}{H^2S^2}\frac{\mathrm{\Lambda }c^2}{3H^2}$$
(6)
Here note that

$$\dot{S}=n\alpha t^\beta ;\beta =n1>0;H=\frac{\dot{S}}{S}=n/t$$
(7)
Using the foregoing relationships, we find from Eq. (6) that for the present epoch $`t=t_1`$ and for the past epoch $`t=t_2`$

$$\frac{\mathrm{\Omega }_11}{\mathrm{\Omega }_21}=\frac{k\alpha ^2t_1^{2\beta }\frac{\mathrm{\Lambda }_1t_1^2c^2}{3}}{k\alpha ^2t_2^{2\beta }\frac{\mathrm{\Lambda }_2t_2^2c^2}{3}}$$
(8)
As $`t_20`$, the second term in the denominator of the above equation $`0`$, and we obtain

$$\frac{\mathrm{\Omega }_11}{\mathrm{\Omega }_21}=\left(\frac{t_1}{t_2}\right)^{2\beta }\left(\frac{\mathrm{\Lambda }_1}{\mathrm{\Lambda }_2}\right)t_1^2t_2^{2\beta }$$
(9)
Again as $`t_20`$, both the terms on the RHS of the foregoing equation $`0`$ if $`\mathrm{\Lambda }`$ varies less sharply than $`t^2`$. Therefore, we find that

$$\mathrm{\Omega }_11=𝒪(0)(\mathrm{\Omega }_21)$$
(10)
So in the present universe, we should not only have $`\mathrm{\Omega }_11𝒪(1)`$, but, we should actually have $`\mathrm{\Omega }_1=\mathrm{\Omega }_{matter}+\mathrm{\Omega }_\mathrm{\Lambda }1`$. Then there is no flatness problem, and on the other hand, the universe should be completely flat at any appreciable value of $`t`$. In fact very recent observations indeed show that this is precisely the case.

## IV Discussion and Conclusion

The age of the universe in the present case could be finite. This can be seen in the following simple way:

$$t=\frac{dS}{\dot{S}}$$
(11)
By using Eq. (7), we see that

$$t=\frac{S^{1/n}}{\alpha }$$
(12)
In contrast there could be more involved models incorporating AE for which one may have $`t=\mathrm{}`$ for $`S0`$. In other words there could be singularity free models for which horizon and flatness problems naturally do not arise. However, we found that in a more orthodox model of AE with a finite age of the universe, horizon and flatness problems cease to exist and, thus, there is no need to postulate, in an adhoc manner the occurrence of a brief spell of inflation. It can be shown that even if the horizon and flatness problems disappear with the idea of inflation in the present epoch, in principle, they may reappear, in future if $`k0`$ and the universe were undergoing DE. On the other hand, such problems do no recur for the present simple model incorporating AE. And although, the singularity and a hot early epoch lie in a finite past, in the present simple case of an orthodox model, note that the universe begins with a “whimper” ($`\dot{S}=0`$ infinity) rather than a “bang”($`\dot{S}=\mathrm{}`$).

We have of course, not suggested any physics which may be driving the universe in the AE mode. This a completely different problem in its own right and many authors are addressing to it. We would simply mention here that it is possible to have practically any form of a time varying $`\mathrm{\Lambda }`$ and consequent desired evolution. And if we consider the AE occuring due to a $`\mathrm{\Lambda }`$ term, which is intrinsic to the fabric of spacetime, the presence of such a term may not significantly alter the local thermodynamics. Thus the relationship between $`S`$ or $`z`$ and temperature could be unaltered in such a scheme. In other words, as an example, the epoch of decoupling may still be at $`z1100`$ as in the standard $`\mathrm{\Lambda }`$-free cosmology. However the $`S(t)`$ or $`z`$ and $`t`$ will change. Consequently the epoch of decoupling will lie at different model dependent cosmic time.


text:
# Energy landscape of a Lennard-Jones liquid: Statistics of stationary points

## Abstract

Molecular dynamics simulations are used to generate an ensemble of saddles of the potential energy of a Lennard Jones liquid. Classifying all extrema by their potential energy $`u`$ and number of unstable directions $`k`$, a well defined relation $`k(u)`$ is revealed. The degree of instability of typical stationary points vanishes at a threshold potential energy $`u_{th}`$, which lies above the energy of the lowest glassy minima of the system. The energies of the inherent states, as obtained by the Stillinger-Weber method, approach $`u_{th}`$ at a temperature close to the mode-coupling transition temperature $`T_c`$.

The complex physical behaviour of supercooled liquids and glasses has always stimulated a description of such systems in terms of dynamical evolution upon a very complicated potential energy landscape . In particular, the presence of many inequivalent glassy minima of the potential energy gives rise to a rich pattern of activated dynamics and has therefore attracted most of the attention of the community in the past years. As a consequence, the analysis of the structure of minima has become the main focus of the energy landscape approach to glasses. Based on the pioneering work of Goldstein  and of Stillinger and Weber , a number of authors have investigated the statistical properties of local minima, which are reached by a gradient descent, starting from an equilibrium configuration of the supercooled liquid state at temperature $`T`$ . The energies of the inherent minima sampled after the quench reveal two characteristic temperatures: At low $`T`$ the system gets trapped in a very small number of basins and falls out of equilibrium. Although this transition temperature depends on the cooling rate, it is very close to the dynamical transition $`T_c`$ of mode-coupling theory (MCT) . A second higher temperature marks the onset of nonexponential relaxation, which has been associated with energy landscape dominated dynamics .

A method which does not focus purely on the properties of minima, is the instantaneous normal mode (INM) approach : the spectrum of the eigenvalues of the Hessian matrix is computed and averaged over all the configurations with the Boltzmann distribution at temperature $`T`$. The INM analysis focuses on two points: First, it has been suggested that the barrier heights and hopping rates can be obtained from the INM spectrum . Second, the temperature where the fraction $`f(T)`$ of negative eigenvalues of the INM spectrum goes to zero is interpreted as the point where the number of directions for free diffusion in phase space vanishes, thus giving an estimate of $`T_c`$, below which activation remains as the only mechanism of diffusion . The INM analysis has been criticized , because equilibrium configurations with unstable directions are in general not stationary points of the potential energy, even if the force vanishes along the unstable directions .

In contrast to the INM method, which is an intrinsically thermal approach by sampling configurations according to their Boltzmann weight, we focus here on the purely geometric properties of the energy landscape. We investigate all the stationary points of the potential energy, be they minima, or unstable saddles . We classify them according to their number of unstable directions, (index K), their energy and the smallest eigenvalue of the Hessian matrix. Thereby we can address the following questions: Is there a threshold energy for saddles, such that for energies below this threshold, it is very unlikely to find saddles and the dynamic relaxation of the system requires activation ? Is there a signature of the threshold energy in the dynamical behaviour of the system ? What is the typical energy difference between saddles of index $`K`$ and $`K+1`$ and can this be taken as an estimate of the potential energy barriers ?

The system under consideration is a binary mixture of large (L) and small (S) particles with $`80\%`$ large and $`20\%`$ small particles. Small and large particles only differ in diameter, but have the same mass. They interact via a Lennard–Jones potential of the form $`V_{\alpha \beta }(r_{ij})=4ϵ_{\alpha \beta }[(\sigma _{\alpha \beta }/r_{ij})^{12}(\sigma _{\alpha \beta }/r_{ij})^6]`$, where $`\stackrel{}{r}_i`$ denotes the position of particle i $`(i=1,2\mathrm{}N)`$ and $`r_{ij}=|\stackrel{}{r}_i\stackrel{}{r}_j|`$. All results are given in reduced units, where $`\sigma _{LL}`$ was used as the length unit and $`ϵ_{LL}`$ as the energy unit. The other values of $`ϵ`$ and $`\sigma `$ were chosen as follows: $`ϵ_{LS}=1.5,\sigma _{LS}=0.8,ϵ_{SS}=0.5,\sigma _{SS}=0.88`$ . The systems were kept at a fixed density $`\rho 1.2`$. Periodic boundary conditions have been applied and the potential has been truncated appropriately according to the minimum image rule . We have applied a truncation procedure which on the one hand ensures the potential energy to be zero at the cut-off $`r_{\mathrm{cut}}`$ and on the other hand provides a continuous first-derivative of $`V_{\alpha \beta }(r)`$ at $`r_{\mathrm{cut}}`$. Throughout this study we present results for systems with $`N=60`$ particles using $`r_{\mathrm{cut}}=1.8`$. A few samples with up to $`N=120`$ have been simulated to validate our results.

In order to explore the geometric properties of the stationary points of the potential energy, we use the following method: We start by equilibrating a random initial configuration of $`N`$ particles at a given temperature $`T`$ using a standard molecular dynamics simulation technique. After equilibration the system evolves for a time $`\tau _{\mathrm{run}}`$. To locate a saddle of the potential energy close to the equilibrated configuration we look for the absolute minima of the modulus square of the force. In order achieve this we perform a quench on a pseudo-potential energy landscape $`W(x)`$ given by $`W(x)=\stackrel{}{}U(x)\stackrel{}{}U(x)`$, whereas the original potential energy is defined by $`U(x)=_{1i<jN}V_{\alpha \beta }(r_{ij})`$. Notice that all absolute minima of $`W(x)`$ are stationary points of $`U(x)`$, hence every saddle of $`U(x)`$ has a well defined basin of attraction. The local minima of $`W(x)`$ do not correspond to zeros of the real force. These points are frequently sampled; they can easily be distinguished from the absolute minima and are excluded from the following analysis.

Given a stationary point we consider the number of negative eigenvalues of its Hessian matrix, that is the index $`K`$ of the saddle. We compute the index density $`k=K/(3N)`$ and the potential energy density $`u=U/N`$ of the saddle and plot these values in the $`(u,k)`$ plane. In Fig.1 we show the results for many saddles, which were sampled with the steepest descent procedure on the pseudo-potential $`W(x)`$ for two different values of the temperature. This plot clearly suggests that there is an underlying curve $`k(u)`$, which is independent of temperature and which encodes a purely geometric feature of the landscape. By sampling stationary points at different values of $`T`$ we are exploring different regions of the potential energy surface and thus different portions of the same geometric curve $`k(u)`$.

To further support this conjecture we have sampled stationary points at various temperatures and have averaged all data in two different ways: Firstly, we have considered all the stationary points with a given index density and we have computed their average potential energy (this is possible because for a finite system $`k`$ can only assume discrete values). We call this procedure geometric average. Secondly, we consider all the stationary points obtained at a given temperature $`T`$ and we compute their average index $`k(T)`$ and energy $`u(T)`$. Eventually, we plot $`k`$ and $`u`$ parametrically in $`T`$ on the $`(u,k)`$ plane. We call this the parametric average. If our sampling of the saddles is a fair exploration of the underlying geometric space, then the two averages must coincide. This is what happens, as it is shown in Fig.2.

The averaged data as plotted in Fig.2 reveal a unique function $`k(u)`$. In other words, if we cut the potential energy landscape with a plane of constant energy density $`u=u_0`$, the stationary points on this plane (or within a narrow shell around this plane) will be dominated by saddles with index density $`k(u_0)`$. Furthermore, $`k(u)`$ is to a very good approximation linear up to an index density of $`10\%`$ negative eigenvalues. This implies that the curve extrapolates to zero at a well defined energy, which we call the threshold energy, $`u_{th}`$, in analogy with spin-glasses . In Fig.3 we present a magnification of the last four points of $`k(u)`$, showing that the linear interpolation of all the data and the linear interpolation of the last four points give the same estimate for the threshold energy, that is $`u_{th}=4.55`$.

The threshold energy marks the border between the saddles-dominated portion of the energy landscape and the minima-dominated one. An interesting point is that $`u_{th}`$ is above the energy of the lowest lying minima we find, $`u_0=4.65`$, as shown in Fig.3 ($`u_0`$ is obtained from an extensive search for minima of the potential energy using the Stillinger-Weber method). This implies the existence of a finite energy density interval where minima are overwhelmingly more numerous than saddles. The same phenomenon has been observed in mean-field models of spin-glasses  exhibiting one step replica symmetry breaking (1RSB).

In Fig.4 we plot the lowest eigenvalue $`\lambda _0`$ of the Hessian as a function of the energy density. As expected, $`\lambda _00`$ for $`uu_{th}`$, implying that most extrema have a small number of eigenvalues close to zero. Somewhat surprising is the approximately linear dependence of $`\lambda _0(u)`$ on energy, $`\lambda _0(u_{th}u)`$, exactly as in 1RSB spin-glasses .

It is important to relate the geometric features of the energy landscape to the dynamical behaviour of the supercooled liquid close to the glass transition. The index density vanishes at the threshold so that minima are the dominant stationary points below $`u_{th}`$. We therefore expect to find a signature of $`u_{th}`$ in the dynamics or, more specifically, a link between the threshold energy and the onset of activated dynamics upon cooling. To that end we have used our MD simulations at temperature T  to compute the energy density $`u_{min}(T)`$ of local minima reached by gradient descent on the potential energy surface. These data are compared in Fig.5 to the difference $`\delta (T)=<U/N>(T)3T/2`$ of the average potential energy density $`<U/N>(T)`$, as calculated in the MD simulation, and the vibrational energy in the harmonic approximation. For a harmonic potential $`\delta (T)`$ is just the energy of the minimum of the well. In Fig.5 we see that $`\delta (T)u_{min}(T)`$ for $`T1.2`$ . This is the range of temperatures which is dominated by the energy landscape and the timescales for the two processes of relaxation - vibrations inside a minimum and hopping between different minima - start to separate. Close to the glass transition the system falls out of equilibrium, as indicated in Fig.5 by the saturation of both quantities, $`u_{min}(T)`$ and $`\delta (T)`$. The temperature where this happens is known to depend on the equilibration time of the MD simulation (see e.g. ).

Extrapolating that part of the curve corresponding to temperatures where equilibration is ensured (see Fig. 5), one observes that both $`u_{min}(T)`$ and $`\delta (T)`$ reach the threshold energy at a temperature approximately equal to the MCT transition $`T_c0.44`$. The crossover from a non-activated dynamics above $`T_c`$ to an activated one below $`T_c`$ , may thus be interpreted in this context as a geometric transition from saddle dominated, to minima dominated regions of the potential energy landscape.

We can obtain an estimate of the energy barrier $`\mathrm{\Delta }U=1/(3k^{})5.0`$ from the slope $`k^{}(u)`$ of the linear function . This estimate has the right order of magnitude, see e.g. . Note however that we do not know which saddles are accessible from a particular minimum via a dynamic trajectory. Furthermore we have not yet computed the entropic contribution to the transition rate, which is necessary for an estimate of the free energy barrier and hence for a comparison with the experimental time scales.

Finally in Fig.6 , we compare our result for $`k(u)`$ to an INM analysis. For the latter we have computed the average fraction of negative eigenvalues for equilibrium configurations as a function of their energy (parametric average). We see that for all temperatures the INM index is higher than the index of genuine saddles at the corresponding energy. Hence the average curvature is overestimated by INM, as expected .

In this Letter we have presented a new method to characterize in a quantitative way the geometric structure of the potential energy landscape, which goes beyond the analysis of local minima. We have shown that the typical saddles index $`k`$ depends linearly on the potential energy density and we have defined a threshold energy, $`u_{th}`$, where the degree of instability vanishes. Furthermore, the threshold energy has been related to the MCT transition temperature $`T_c`$.

The method of steepest descent on the pseudo-potential $`W(x)`$ is powerful enough to provide detailed information on the properties of saddles, so that one can estimate not only energy barriers but also the entropic contribution to the transition rates. Given that we know all eigenvalues of the Hessian at the stationary points, we can evaluate the partition function at the saddle, which is needed as an input for the transition state theory. So far we have mainly looked at one system size, i.e. $`60`$ particles, and have only checked a few results for larger systems. It would be interesting to study the effects of varying system size on barrier energies, entropies and number of unstable directions. Another important point is the mutual accessibility of minima and saddles. This latter problem may be tractable by a combination of both, steepest descent on the potential energy and steepest descent on the pseudo-potential surface. Work along these lines is in progress.

###### Acknowledgements.

A.C. thanks A. Bray, F. Colaiori, T. Keyes, M. Moore, G. Parisi, F. Ritort, S. Sastry and F. Thalmann for useful discussions. The work of A.C. was supported by EPSRC(UK) Grant GR/L97698 and the work of K.B. by the DFG under grant No. Br 1894/1-1. In memory of Kurt.


text:
# Bell State Preparation using Pulsed Non-Degenerate Two-Photon Entanglement

## Abstract

We report a novel Bell state preparation experiment. High-purity Bell states are prepared by using femtosecond pulse pumped *nondegenerate* collinear spontaneous parametric down-conversion. The use of femtosecond pump pulse does not result in reduction of quantum interference visibility in our scheme in which post-selection of amplitudes and other traditional mechanisms, such as, using thin nonlinear crystals or narrow-band spectral filters are not used. Another distinct feature of this scheme is that the pump, the signal, and the idler wavelengths are all distinguishable, which is very useful for quantum communications.

Preparation and measurement of the Bell states are two important issues in modern quantum optics, especially for quantum communications, quantum teleportation, etc . For photons, such states can be realized by using the entangled photon pairs generated in spontaneous parametric down-conversion (SPDC). By making appropriate local operations on the SPDC photon pairs, one can prepare all four Bell states.

The polarization Bell states, for photons, can be written as

$`|\mathrm{\Phi }^\pm `$ $`=`$ $`|X_1,X_2\pm |Y_1,Y_2,`$ (1)
$`|\mathrm{\Psi }^\pm `$ $`=`$ $`|X_1,Y_2\pm |Y_1,X_2,`$ (2)
where the subscripts 1 and 2 refer to two different photons, photon 1 and photon 2, respectively, and they can be arbitrarily far apart from each other. $`|X`$ and $`|Y`$ form the orthogonal basis for the polarization states of a photon, for example, it can be horizontal ( $`|H`$) and vertical ($`|V`$) polarization state, as well as $`|45^{}`$ and $`|45^{}`$, respectively. This means that the quantum interference should be independent of the choice of the bases.

Such an experiment was first performed by Shih and Alley in which non-collinear type-I SPDC and a beamsplitter were used to prepare a Bell state , but it is very difficult to align such a system. Collinear type-II SPDC is thus developed . There is, however, a common problem: the entangled photon pairs have 50% chances of leaving at the same output ports of the beamsplitter. Therefore, the state prepared after the beamsplitter may not be considered as a Bell state without amplitude post-selection . Only when one considers the coincidence contributing terms by throwing away two out of four amplitudes (post-selection of 50% of the amplitudes), the state is then said to be a Bell state. This problem is later solved by using non-collinear type-II SPDC or using two non-collinear type-I SPDC .

In the cw pumped SPDC, entangled photon pairs occur randomly since the process is “spontaneous”, so whereabouts of the photon pair is completely uncertain within the coherence length of the pump laser beam. This huge time uncertainty makes it difficult for applications such as generation of multi-photon entangled state, quantum teleportation, etc, as interactions between entangled photon pairs generated from different sources are required. This difficulty was thought to be solved by using a femtosecond pulse laser as a pump. Unfortunately, femtosecond pulse pumped type-II SPDC shows poor quantum interference visibility due to the very different (compared to the cw case) behavior of the two-photon effective wave-function . One has to utilize special experimental schemes to achieve complete overlap of the two-photon amplitudes. Traditionally, the following methods were used to restore the quantum interference visibility in femtosecond pulse pumped type-II SPDC: (i) use a thin nonlinear crystal ($`100\mu `$m)  or (ii) use narrow-band spectral filters in front of detectors . Both methods, however, reduce the available flux of the entangled photon pair significantly  and cannot achieve complete overlap of the wave-functions in principle .

The first attempt to achieve high-visibility quantum interference in femtosecond pulse pumped type-II SPDC without using narrowband filters and a thin crystal was reported in Ref.. The observed visibility, however, was rather low and keeping the phase coherence over a long term would be very difficult since a Michelson interferometer is used. Also, such a scheme cannot be used to prepare a Bell state. Recently, we reported a high-visibility quantum interference experiment in which photon pairs are entangled both in polarization and space-time using femtosecond pulse pumped type-I SPDC . However, it cannot be considered as a true Bell state preparation since post-selecting 50% of the amplitudes was still necessary.

In this Letter, we report a Bell state preparation experiment in which we effectively eliminate *any post-selection* in femtosecond pulse pumped SPDC for the first time. Other features in our scheme include: (i) collinear SPDC makes the alignment much easier, (ii) Alice and Bob shares photons of different frequencies entangled in both space-time and polarization, (iii) phase coherence is automatically kept and the visibility as high as 92% is observed, (iv) thick crystals can be used to increase the intensity (without losing the visibility), and (v) the spectral bandwidth is reduced significantly by the use of *nondegenerate* SPDC. These features make our scheme a good source of Bell states for quantum information experiments.

The basic idea of the experiment is illustrated in Fig.1(a). A $`45^{}`$ polarized femtosecond laser pulse (central wavelength $`\lambda _p=400`$nm and pulse duration $`\sigma =80`$fsec.) enters the Mach-Zehnder interferometer (MZI) which contains a type-I nonlinear crystal in each arm. One crystal has its optic axis oriented vertically ($``$) and another horizontally ($``$). Polarizing beamsplitter (PBS) splits the $`45^{}`$ polarized pump pulse into the vertical and horizontal polarized pulses propagating along different arms of the MZI. Then non-degenerate collinear type-I SPDC occurs, with equal probability, at each crystal (signal wavelength $`=730`$nm and idler wavelength $`=885`$nm) and they are mixed at the dichroic beamsplitter which directs 730nm photons to detector $`D_1`$ and 885nm to detector $`D_2`$. In the simplified single mode approximation, the quantum state generated from the vertically oriented crystal ($``$) is $`|\psi _1=|H_{730}|H_{885}`$ and from the horizontally oriented one ($``$) is $`|\psi _2=|V_{730}|V_{885}`$. $`H`$ and $`V`$ represent horizontal and vertical polarization state of a single photon respectively. Subscripts 730 and 885 refer to the wavelengths 730nm and 885nm, respectively. When the MZI is balanced, the quantum state after the MZI is (without throwing away any amplitudes)

$$|\mathrm{\Phi }=|V_{730}_1|V_{885}_2+e^{i\mathrm{\Delta }\phi }|H_{730}_1|H_{885}_2,$$
(3)
where $`\mathrm{\Delta }\phi `$ is the relative phase between the two amplitudes and it can easily be varied by scanning one of the mirror of the MZI.

The coincidence counting rate ($`R_c`$) is calculated as ,

$$R_c=𝑑t_+𝑑t_{12}\left|𝒜(t_+,t_{12})\right|^2,$$
(4)
where $`t_{12}=t_1t_2`$ and $`t_+=(t_1+t_2)/2`$. $`t_i=T_il_i/c`$ where $`T_i`$ is the time at which detector $`i`$ fires and $`l_i`$ is the optical path length from the surface of the crystal to the detector $`i`$. $`𝒜(t_+,t_{12})`$ is the amplitude of the biphoton as explicitly calculated in Ref..

For the scheme shown in Fig.1(a), $`𝒜(t_+,t_{12})`$ is the sum of the two amplitudes originated from the crystal in each arm of the MZI:

$$𝒜(t_+,t_{12})=𝒜_a(t_+,t_{12})+𝒜_b(t_+,t_{12}),$$
(5)
where the subscripts $`a`$ and $`b`$ refer to the crystal from which the amplitudes are created. The delay $`T`$ introduced in one arm modifies the amplitude $`𝒜_b(t_+,t_{12})𝒜_b(t_++T,t_{12})`$ and determines the additional phase shift for the biphoton amplitudes $`\mathrm{\Delta }\phi =\mathrm{\Omega }_pT=K_p\mathrm{\Delta }x`$, where $`K_p=2\pi /\mathrm{\Lambda }_p`$, $`\mathrm{\Omega }_p(\mathrm{\Lambda }_p)`$ the central frequency (wavelength) of the pump, and $`\mathrm{\Delta }x`$ the spatial delay. Due to the energy conservation and negligibly small dispersion of the air, the phase shift depends only on the pump wavelength $`\mathrm{\Delta }\phi =K_s\mathrm{\Delta }x+K_i\mathrm{\Delta }x=K_p\mathrm{\Delta }x=\mathrm{\Delta }\phi _p`$ although the delay is introduced to the SPDC field . If the crystals are the same and the pump fields in different arms of the MZI are identical,

$$|𝒜_a(t_+,t_{12})|=|𝒜_b(t_+,t_{12})|.$$
(6)
The coincidence counting rate is then calculated to be

$$R_c=1+V\mathrm{cos}(\mathrm{\Omega }_pT),$$
(7)
where $`V1`$ in this experiment . Note that the angles of the analyzers $`A_1`$ and $`A_2`$ are assumed to be $`45^{}`$. From Eq.(7), we expect that the coincidence counting rate will be modulated in the pump central wavelength when $`T`$ is varied.

There are also two more ways to vary the phases of interference by introducing relative delays (using a piece of birefringent material, such as a quartz plate) after the output beamsplitter, i.e., in the signal ($`\mathrm{\Delta }\phi _s`$) and/or in the idler ($`\mathrm{\Delta }\phi _i`$) channels. Therefore we obtain

$$R_c=1+V\mathrm{cos}(\mathrm{\Delta }\phi _p\mathrm{\Delta }\phi _i\mathrm{\Delta }\phi _s),$$
(8)
where $`\mathrm{\Delta }\phi _p`$, $`\mathrm{\Delta }\phi _i`$, and $`\mathrm{\Delta }\phi _s`$ refer to the relative phases of the pump, the idler, and the signal, respectively.

As we have shown so far, one can eliminate the possibility of the entangled photon pairs leaving at the same output ports of the beamsplitter by employing non-degenerate two-photon entanglement. In this scheme, high-visibility quantum interference can be achieved independent of the crystal thickness and the spectral filter bandwidths even with femtosecond pulse pump.

In practice, however, one would not like to use a MZI in the experimental setup due to stability related issues. Therefore, we use a collinear scheme where two type-I BBO crystals are placed collinearly in the pump beam path, see Fig1(b). Two type-I BBO crystals with thickness 3.4mm each (the first one is oriented horizontally and the second one is oriented vertically) are then pumped by a $`45^{}`$ polarized pump pulse. As described before, the quantum state resulting from the first BBO is $`|V_{730}|V_{885}`$ and that from the second BBO is $`|H_{730}|H_{885}`$. Since both crystals are pumped equally, the two amplitudes are equally probable. Due to the dispersion, however, $`|V_{730}|V_{885}`$ from the first BBO ($``$) and $`|H_{730}|H_{885}`$ from the second BBO ($``$) are distinguishable in time. To make $`|V_{730}|V_{885}`$ and $`|H_{730}|H_{885}`$ indistinguishable in time, one needs to compensate the delay experienced by the SPDC photon pairs at each crystal. This compensation can be made by using a properly oriented quartz rod. If the compensation is made properly, either before or after the down-conversion non-linear crystals, one will observe high-visibility quantum interference without any spectral post-selection.

In the collinear scheme, having a perfect temporal compensation is difficult when the signal wavelength differs very much from the idler wavelength. This is because the signal-idler photon pairs created from the first BBO ($``$) experience different dispersion when they pass through the second BBO ($``$). (The MZI scheme does not have this disadvantage). In this experiment, for the wavelengths we are interested in, the temporal separation is rather small and it does not affect the interference visibility. To prevent further dispersion effects, the compensation is made before the BBO crystals. The compensator consists of a quartz rod and two quartz plates whose optic axes are oriented vertically, see Ref., and it imposes roughly 1.5psec required delay between the $`H`$\- and $`V`$-polarized 400nm pump pulse which is mainly determined by the thickness of the BBO crystals. By tilting the two quartz plates in the opposite directions, the phase delay $`\mathrm{\Delta }\phi _p`$ can be varied to prepare a Bell state. After the two BBO crystals, the remaining UV radiation is blocked by a UV reflecting mirror and the collinear SPDC is selected by a diaphragm. Then a dichroic beamsplitter is used to reflect the signal (730nm) to $`D_1`$ and to transmit the idler (885nm) to $`D_2`$. Two quartz plates are inserted in each beam path to vary the relative phase of the signal or the idler independently. The detector package consists of a single-photon counting module, an interference filter which is used to cut the pump noise , and a polarization analyzer.

To demonstrate the effectiveness of this scheme, we first study the space-time interference as a function of $`\mathrm{\Delta }\phi _p`$ by setting $`\theta _1=\theta _2=45^{}`$, where $`\theta _1`$ and $`\theta _2`$ are the angles of the analyzers $`A_1`$ and $`A_2`$ (measured from the vertical direction). According to Eq.(8), one should observe a pump wavelength modulation in the coincidence counting rate. Note that $`\mathrm{\Delta }\phi _s`$ and $`\mathrm{\Delta }\phi _i`$ are fixed. The observed modulation period is 400nm, see Fig.2(a), which agrees with the theory.

To prepare $`\mathrm{\Phi }^+`$ ($`\mathrm{\Phi }^{}`$) state, identified by constructive (destructive) interference, one just needs to set $`\mathrm{\Delta }\phi _p\mathrm{\Delta }\phi _s\mathrm{\Delta }\phi _i=0,2\pi ,4\pi ,\mathrm{}`$ ($`\mathrm{\Delta }\phi _p\mathrm{\Delta }\phi _s\mathrm{\Delta }\phi _i=\pi ,3\pi ,5\pi ,\mathrm{}`$) which can be done by tilting the quartz plates so that the space-time interference fringe is at the maximum (minimum). Note that $`\mathrm{\Psi }^+`$ and $`\mathrm{\Psi }^{}`$ Bell states can also be easily prepared by introducing a $`\lambda /2`$ plate in one output port of the dichroic beamsplitter.

We have also experimentally demonstrated the polarization interference for $`\mathrm{\Phi }^+`$ and $`\mathrm{\Phi }^{}`$. For $`\mathrm{\Phi }^\pm `$, the coincidence counting rate is calculated to be

$$R_c\left|\theta _2,\theta _1|\mathrm{\Phi }^\pm \right|^2\mathrm{cos}^2(\theta _1\theta _2).$$
(9)
This means that one should observe high-visibility modulation in polarization correlation measurement for arbitrary values of $`\theta _1`$ and $`\theta _2`$. To confirm this experimentally, we first set $`\theta _1=45^{}`$ and varied $`\theta _2`$. High-visibility polarization correlation is observed, see Fig.2(b). We then repeated this measurement for many other values of $`\theta _1`$ and observed that the visibility remained the same. This confirms Eq.(9). In other words, we have successfully prepared polarization Bell states from femtosecond pulse-pumped SPDC without amplitude and spectral post-selection.

Unlike the usual degenerate two-photon sources, this source has one distinctive feature: one can vary three different phases independently which is very useful for quantum communications. To demonstrate this interesting feature, we observe the space-time interference by varying the relative phases of the signal ($`\mathrm{\Delta }\phi _s`$) and the idler ($`\mathrm{\Delta }\phi _i`$) independently. In these measurements, $`\theta _1`$ and $`\theta _2`$ are set at $`45^{}`$. The effect of $`\mathrm{\Delta }\phi _p`$ is already demonstrated in Fig.2(a). When the signal phase $`\mathrm{\Delta }\phi _s`$ is varied, see Fig.3(a), the signal wavelength (730nm) modulation is observed in coincidence rate, while varying the idler phase $`\mathrm{\Delta }\phi _i`$, see Fig.3(b), the idler wavelength (885nm) modulation is observed. Fig. 2(a), Fig.3(a), and Fig.3(b) clearly demonstrate Eq.(8).

Finally, we tested the condition $`\mathrm{\Delta }\phi _p=\mathrm{\Delta }\phi _s+\mathrm{\Delta }\phi _i`$ by varying the signal phase $`\mathrm{\Delta }\phi _s`$ and the idler phase $`\mathrm{\Delta }\phi _i`$ at the same time. The quartz plates in both the signal and the idler paths are tilted at the same time with equal angles. To see whether the data agrees with the theory, total phases accumulated in both beam path are calculated from the tilt angle, i.e., $`\mathrm{\Delta }\phi _s+\mathrm{\Delta }\phi _i`$. As evidenced from Fig.3(c), the data agrees well with Eq.(8).

In summary, we have demonstrated a novel scheme to prepare pulsed entangled photon pairs from which all four Bell states can be easily obtained. Amplitude and spectral post-selection are not necessary. Note also that non-maximally entangled states can be prepared by changing the relative intensities of the pump beams. The visibility and the photon flux is greatly enhanced by this method although a femtosecond pulse laser is used as a pump. The signal, the idler, and the pump phases can be varied independently with different modulation frequency which is very useful for quantum communications.

We would like to thank M.H. Rubin for helpful discussions. This research was supported, in part, by the Office of Naval Research, ARDA and the National Security Agency. S.P.K. also thanks the Russian Fund for Fundamental Research Grant No. 99-02-16419 for partially supporting his visit to UMBC.