text:
# Ferromagnetism in the Hubbard model with long-range and correlated hopping

## Abstract

The extrapolation of small-cluster exact-diagonalization calculations is used to examine ferromagnetism in the one-dimensional Hubbard model with long-range and correlated hopping. It is found that the correlated hopping term stabilizes the ferromagnetic state for a wide range of electron interactions $`U`$ and electron concentrations $`n`$. The critical value of the interaction strength $`U_c`$ above which the ferromagnetic state becomes stable is calculated numerically and the ground-state phase diagram of the model is presented for physically the most interesting cases.

The Hubbard model has become, since its inception  in 1963, one of the most popular examples of a system of interacting electrons with short-range interactions. It has been used in the literature to study a great variety of many-body effects in metals, of which ferromagnetism, metal-insulator transitions, charge-density waves and superconductivity are the most common examples. Of all these cooperative phenomena the problem of ferromagnetism in the Hubbard model has the longest history. Although the model was originally introduced to describe the band ferromagnetism of transition metals, it soon turned out that the single-band Hubbard model is not the canonical model for ferromagnetism. In fact the existence of saturated ferromagnetism has been proven rigorously only for very special limits. The first well-known example is the Nagaoka ferromagnetism that comes from the Hubbard model in the limit of infinite repulsion and one hole in a half-filled band . Another example where saturated ferromagnetism has been shown to exist is the case of the one-dimensional Hubbard model with nearest and next-nearest-neighbor hopping at low electron densities . Moreover, several examples of the fully polarized ground state have been found on special lattices (special conduction bands) as are the fcc-type lattices , the bipartite lattices with sublattices containing a different number of sites , the lattices with unconstrained hopping of electrons  and the flat bands . This indicates that the lattice structure and the kinetic energy of electrons, i.e., the type of hopping play an important role in stabilizing the ferromagnetic state.

In our previous paper  we have shown that if the electron hopping is described by a more realistic model (than the nearest-neighbor hopping) then ferromagnetism comes from the Hubbard model naturally for a wide range of the model parameters. No extra interactions terms should be included. In particular, we have found that the long-range hopping with power decaying hopping amplitudes $`t_{ij}`$ given by

$$t_{ij}(q)=\{\begin{array}{cc}0,\hfill & i=j,\hfill \\ q^{|ij|}/q,\hfill & ij\hfill \end{array}$$
(1)
gives rise to ferromagnetism for electron densities above half-filling. As soon as $`q`$ that controls the effective range of the hopping ($`0q1`$) is different from zero the ferromagnetic state is stabilized for all Coulomb interactions $`U`$ greater than some critical interaction strength $`U_c(q)`$ that value is dramatically reduced with increasing $`q`$. From this point of view one of main reasons why the ferromagnetic state absent in the ordinary Hubbard model with nearest-neighbor hopping ($`q=0`$) is that the description of the electron hopping was too simplified.

Here we further generalize this model by introducing the correlated hopping term, in which the $`\sigma `$-electron hopping amplitudes between lattice sites $`i`$ and $`j`$ depend explicitly on $`n_{i\sigma }`$ and $`n_{j\sigma }`$ occupancy, i.e.

$$t_{ij}^\sigma =t_{ij}[1+t^{}(n_{i\sigma }+n_{j\sigma })].$$
(2)
The importance of the correlated hopping term on the ground-state properties of the Hubbard model has been already mentioned by Hubbard . Later Hirsch  pointed out that this term may be relevant in explanation of superconducting properties of strongly correlated electrons. Here we examine effects of this term on the stability of the fully polarized ferromagnetic state. The same subject has been studied recently by Amadon and Hirch , as well as by Kollar and Vollhardt , however they considered hopping only between the nearest-neighbor sites.

The selection of hopping matrix elements in the form given by Eq. 1 and Eq. 2 has several advantages. It represents a much more realistic type of electron hopping on a lattice (in comparison to nearest-neighbor hopping), and it allows us to change continuously the type of hopping and thus immediately study the effect of long-range and correlated hopping on the ground state properties of the Hubbard model.

The Hamiltonian of the single-band Hubbard model in which the effects of long-range and correlated hopping are incorporated is given by

$$H=\underset{ij\sigma }{}t_{ij}^\sigma c_{i\sigma }^+c_{j\sigma }+U\underset{i}{}n_in_i,$$
(3)
where $`c_{i\sigma }^+`$ and $`c_{i\sigma }`$ are the creation and annihilation operators for an electron of spin $`\sigma `$ at site $`i`$, $`n_{i\sigma }`$ is the corresponding number operator ($`N_\sigma =_in_{i\sigma }`$) and $`U`$ is the on-site Coulomb interaction constant.

To examine possibilities for existence of ferromagnetism in this model the ground states are determined by exact diagonalizations for a wide range of model parameters ($`q,t^{},U,N=_\sigma N_\sigma `$). Typical examples are then chosen from a large number of available results to represent the most interesting cases. The results obtained are presented in the form of phase diagrams in the $`U`$-$`q`$ plane. To determine the phase diagram in the $`U`$-$`q`$ plane (corresponding to some $`t^{}`$ and $`N`$) the ground state energy of the model is calculated point by point as functions of $`q`$ and $`U`$. Of course, such a procedure demands in practice a considerable amount of CPU time, which imposes severe restrictions on the size of clusters ($`L`$) that can be studied with this method ($`L16`$). Fortunately, we have found that the ground-state energy of the model depends on $`L`$ only very weakly (for a wide range of the model parameters) and thus already such small clusters can describe satisfactorily the ground state properties of the model.

The results of our small-cluster exact-diagonalization calculations obtained on finite clusters up to $`L=16`$ sites are summarized in Fig. 1 and Fig. 2. There is shown the critical interaction strength $`U_c`$, above which the ground state is ferromagnetic, as a function of $`q`$ for selected values of $`n=N/L`$ and $`t^{}`$ ($`n=3/2,7/4;t^{}=0,0.2,0.4`$). To reveal the finite-size effects on the stability of ferromagnetic domains, the behavior of the critical interaction strength $`U_c(q)`$ has been calculated on several finite clusters at each electron filling. It is seen that finite-size effects on $`U_c`$ are small and thus these results can be satisfactorily extrapolated to the thermodynamic limit $`L\mathrm{}`$. Our results clearly demonstrate that the ferromagnetic state is strongly influenced by correlated hopping ($`t^{}`$) and generally it is stabilized with increasing $`t^{}`$. The effect is especially strong for intermediate and strong values of $`q`$. Even, there exists some critical value of $`q`$ above which the ground state is ferromagnetic for all nonzero $`U`$. With increasing $`t^{}`$ this critical value shifts to lower values of $`q`$ (that represent a much more realistic type of electron hopping) and the ferromagnetic domain correspondingly increases. Performing exhaustive numerical studies of the model for a wide range of electron concentrations (on different lattice clusters) we have found that the model exhibits the same behavior for all electron concentrations above half-filling , and that with increasing concentration this effect becomes more pronounced (see Fig. 1 and Fig. 2). These results clearly show that ferromagnetism comes naturally from the Hubbard model with long-range and correlated hopping for a wide range of model parameters without any other assumptions. This opens a new route towards the understanding of ferromagnetism in the Hubbard model.

In summary, the extrapolation of small-cluster exact-diagonalization calculations was used to examine ferromagnetism in the one-dimensional Hubbard model with long-range and correlated hopping. It was found that the correlated hopping term stabilizes the ferromagnetic state for a wide range of electron interactions $`U`$ and electron concentrations $`n`$. The critical value of the interaction strength $`U_c`$ above which the ferromagnetic state becomes stable was calculated numerically and the ground-state phase diagram of the model is presented for physically the most interesting cases.

This work was supported by the Slovak Grant Agency VEGA under Grant No. 2/7021/20 and the Science and Technology Assistance Agency under Grant APVT-51-021602. Numerical results were obtained using computational resources of the Computing Centre of the Slovak Academy of Sciences.

Figure Captions

Fig. 1. The critical interaction strength $`U_c`$ as a function of $`q`$ calculated for different $`t^{}`$ and $`L`$ at $`n=3/2`$. Curves from up to down correspond to: $`t^{}=0,0.2`$ and 0.4.

Fig. 2. The critical interaction strength $`U_c`$ as a function of $`q`$ calculated for different $`t^{}`$ and $`L`$ at $`n=7/4`$. Curves from up to down correspond to: $`t^{}=0,0.2`$ and 0.4.


format: text

text:
# Periodic orbit theory of two coupled Tchebyscheff maps

## 1 Introduction

Coupled map lattices (CML) are dynamical systems with discrete space and time and continuous state variables. Their most common form in applications is the diffusively coupled model

$$\mathrm{\Phi }_{n+1}^{(i)}=(1a)f(\mathrm{\Phi }_n^{(i)})+(a/2)[g(\mathrm{\Phi }_n^{(i+1)})+g(\mathrm{\Phi }_n^{(i1)})]$$
(1)
where $`n`$ is a discrete time step, $`i`$ is a lattice point with a periodic boundary condition, $`f(\mathrm{\Phi })`$ the local map, $`g(\mathrm{\Phi })`$ the coupling and $`a`$ a continuous parameter. They prove particularly useful in the investigation of chaotic dynamics of spatially extended systems (spatiotemporal chaos), but have also been employed in chemical, biological and engineering modelling, as well as pattern formation. The discrete structure in space and time can either be the approximation of a less tractable continuous system (as in turbulence), or may arise naturally in applications (as in population dynamics of discrete colonies of organisms). Refs.  give a good overview of early work to 1998 including the above applications; recent developments include advances in the mathematics of the relevant Banach spaces and transfer operators  as well as ever widening applications; examples include ecology , neurons , traffic flow  cryptography  and field theory .

From the point of view of understanding the dynamics of coupled chaotic maps, the most common CML studied are where the local dynamics $`f`$ is quadratic (the “logistic” map); this class of one dimensional maps is now reasonably well understood. The simplest example of a quadratic map is the second degree Tchebyscheff polynomial $`f(\mathrm{\Phi })=T_2(\mathrm{\Phi })=2\mathrm{\Phi }^21`$ which is well known to be conjugated to a piecewise linear map with complete binary symbolic dynamics, and hence exactly solvable.

In this paper we consider diffusively coupled Tchebyscheff maps (see Sec. 2.2 for details) of degree 2 and 3. These have been previously used by Beck  in some interesting but unconfirmed ideas concerning dynamical generation of noise for stochastic quantisation leading to predictions of parameters in the standard model of particle physics. One of the authors  subsequently showed that four of Beck’s six CML models exhibit stable synchronised time periodic states at arbitrarily small values of the coupling parameter $`a`$, albeit with small basins of attraction. The motivation for considering these models derives not only from possible applications, but also theoretical simplicity in that the uncoupled state is fully chaotic and exactly solvable, the small coupling limit has interesting structure, and the equations (being low degree polynomials) allow many exact calculations.

Here we apply classical periodic orbit theory  to a lattice of two coupled Tchebyscheff maps. This theory allows very precise computations of long time properties (averages, Lyapunov exponents, dimensions, decay of correlations, etc.) of chaotic systems, but has been mostly restricted to low dimensional systems in the past. This paper contains the first steps towards developing a high dimensional periodic orbit theory, making use of the properties of Tchebyscheff maps to locate a large number of periodic orbits. If the Beck theory turns out to be true, the insight gained regarding these CML will be invaluable in developing periodic orbit based computation methods which are potentially far more precise than direct simulation, and which are required in particle physics. More generally, a periodic orbit theory of CML will allow a far greater precision and insight into the many spatiotemporal systems discussed at the beginning than is currently available using direct methods.

In order to achieve our immediate goal of computing Beck’s observables on a two site lattice, we first need to look for periodic orbits of (4) by using appropriate numerical methods and study some of their properties, as described in the following sections. The locations, the stability and the number of periodic orbits of such maps depend on the value of $`a`$ and current results do not allow prediction of these changes in detail. As we discuss in section 2, unexpected phenomena take place as soon as the coupling parameter $`a`$ is increased from 0, for example a few periodic points suddenly move off the region $`\mathrm{\Phi }^i[1,1]`$, which we consider for our analysis. As $`a`$ is increased from, periodic points drift on the domain $`\mathrm{\Phi }^i[1,1]`$, moving into the complex plane and then reappearing on the real axis as $`a`$ is increased further. The stability of periodic orbits is also heavily affected by the value of the coupling and stable and/or almost stable orbits appear for certain $`a`$’s. Such behaviour has important consequences for the convergence of cycle expansions , as discussed in section 3. In fact both the rate and the kind of convergence of any average we compute on the lattice strongly depends on the value of the coupling parameter. Still in the same section we address the applications of equations (4) to the theory of Beck. As pointed out by one of the authors , the system (4) does not show ergodicity for all values of $`a`$, therefore direct simulation is not always guaranteed to be reliable in computing averages over the lattice. That is one reason why we investigate the efficiency of averaging over the periodic points. When it is possible to achieve a significant precision for the zeta function (see sec. 3.1), we have computed the average of some string potentials defined in sec. 3.2, and then compared the results with direct simulation on the same lattice. Conclusions and outlook are presented in section 4.

## 2 Periodic solutions

### 2.1 Finding periodic orbits

The first problem is to find periodic solutions to the map (4) up to the highest possible period. A number of numerical methods can be employed for that purpose. In our case the quickest one seems to be the inverse iteration method, as follows: Let $`x^{}=f(x)`$ be a map whose cycles are all expanding111they have no marginal or contracting directions, then it follows that its inverse $`x^{}=f^1(x)`$ has only stable periodic orbits, which can easily be found by iteration. A different symbol sequence (where symbols are associated with branches of the multivalued inverse) is to be used for every orbit . Unfortunately this method cannot be applied a priori, unless we are certain that the map has only expanding cycles. It turns out that the system (4) has got a few periodic orbits whose Jacobian has one stable and one unstable eigenvalue (see sec 2.3). This makes it impossible for inverse iteration to find them, since we cannot guess initial conditions that lie right on the unstable manifold.

A different approach to the problem is that to find the zeros of the equation $`f^n(x)x=0`$, where $`n`$ is the period in question. The multipoint Newton’s method  is believed to be one of the best algorithms for that purpose. It has the advantage of converging relatively fast to the solution, provided that sufficiently accurate initial conditions are given. A cycle of period $`n`$ is a root of the vector function $`F`$:

$$F(x)=F\left(\begin{array}{c}x_1\\ x_2\\ \mathrm{}\\ x_n\end{array}\right)=\left(\begin{array}{c}x_1f(x_n)\\ x_2f(x_1)\\ \mathrm{}\\ x_nf(x_{n1})\end{array}\right)$$

The iteration we need in order to approach the roots has the form

$$\frac{d}{dx}F(x)(x^{}x)=F(x)x^{}=x\left[\frac{d}{dx}F(x)\right]^1F(x)$$

Therefore, after a certain number of iterations we expect to find one root. Of course the method has to be performed every time with different initial conditions, until we are sure all the solutions have been found. The issue concerning the choice of initial conditions is a delicate one, because in principle it is difficult to guess the initial points on a 2-dimensional plane accurately. Nevertheless, it is useful to remember that we are dealing with weakly coupled maps ($`a1`$). For that reason we expect the periodic points to be relatively close to those of the uncoupled system ($`a=0`$), that can be found analytically (see sec. 2.2) and hence exactly.

### 2.2 Uncoupled system

The spatially uncoupled dynamics of Tchebyscheff maps is given by

$$\mathrm{\Phi }_{n+1}^i=T_N(\mathrm{\Phi }_n^i)$$
(2)
where $`T_N`$ is the $`N`$th Tchebyscheff polynomial, defined in general as

$$T_N(\mathrm{\Phi })=\mathrm{cos}(N\mathrm{arccos}\mathrm{\Phi })$$
(3)
of which we consider the first two, namely

$$T_2(\mathrm{\Phi })=2\mathrm{\Phi }^21$$

$$T_3(\mathrm{\Phi })=4\mathrm{\Phi }^33\mathrm{\Phi }$$

The general (coupled) Tchebyscheff maps have the form

$$\mathrm{\Phi }_{n+1}^i=(1a)T_N(\mathrm{\Phi }_n^i)+s\frac{a}{2}(T_N^b(\mathrm{\Phi }_n^{i1})+T_N^b(\mathrm{\Phi }_n^{i+1}))$$
(4)
where $`a`$ is the string coupling parameter, $`s`$ can be $`\pm 1`$ and $`b`$ can take on the values 1 and 0, determining respectively advanced coupling ($`T_N^1(\mathrm{\Phi })=T_N(\mathrm{\Phi })`$, “A”) or backward coupling ($`T_N^0(\mathrm{\Phi })=\mathrm{\Phi }`$, “B”). We have then six models $`2A^+`$, $`2A^{}`$, $`2B^+`$, $`2B^{}`$, $`3A^+`$ and $`3B^+`$. The other two $`N=3`$ maps are equivalent to these by a trivial symmetry.

The model (4), becomes

$`x_{n+1}=(1a)T_N(x_n)+saT_N^b(y_n)`$ (5)
$`y_{n+1}=(1a)T_N(y_n)+saT_N^b(x_n)`$ (6)
on a two site lattice. In this section we consider the uncoupled case $`a=0`$, in which the system (5,6) reduces to

$$x^{}=T_N(x)y^{}=T_N(y)$$
(7)
This is useful as there is a close correspondence between orbits at $`a=0`$ and orbits with $`a>0`$. This correspondence is used to find orbits in the coupled dynamics using the uncoupled orbits as starting points for the multipoint Newton algorithm. The correspondence is not exact — the $`N^{2n_p}`$ degree polynomial equation giving the periodic condition always has this number of (possibly degenerate) solutions, however for $`a>0`$ some of these solutions are not in the range $`[1,1]`$ or are complex, and hence not relevant to the dynamics. This is discussed further in the next section.

We begin first with some definitions. A periodic point of length $`n_p`$ is a solution of $`f^{n_p}(x)=x`$; as discussed above, there are at most $`N^{2n_p}`$ of them. A periodic orbit is the orbit $`\{f^j(x)\},j=0,1,\mathrm{}n_p1`$ of a periodic point. A prime orbit is a periodic orbit which is not the repeat of a shorter orbit, that is, $`f^j(x)=x`$ only if $`j`$ is a multiple of $`n_p`$. $`(x,y)`$ is a periodic point of length $`n_p`$ of the uncoupled system if and only if $`x`$ and $`y`$ are both periodic points of length $`n_p`$ of the map $`T_N`$. It is well known  that the $`T_N`$ maps are conjugate to piecewise linear maps with complete $`N`$-ary symbolic dynamics, by the transformation $`x=\mathrm{cos}\pi \theta `$. This leads to exact expressions for all the periodic points of length $`n_p`$, namely $`x=\mathrm{cos}[2j\pi /(N^{n_p}\pm 1)]`$ with $`0jN^{n_p}/2`$. From these we enumerate all the prime orbits of the uncoupled two dimensional dynamics, eliminating cyclic permutations which correspond to points on the same orbit, and repeats of shorter orbits.

As an example we consider $`N=n_p=2`$. Let symbol 0 (1) denote $`x[1,0]`$ (respectively $`x(0,1]`$); in general there are $`N`$ symbols. In the symbolic representation of periodic points, the spatial variable ($`i`$ in Eq. (1)) will increase to the right, and the temporal variable ($`n`$ in Eq. (1)) will increase downwards. The periodic points of the one dimensional map are

$$\begin{array}{c}0\\ 0\end{array},\begin{array}{c}0\\ 1\end{array},\begin{array}{c}1\\ 0\end{array},\begin{array}{c}1\\ 1\end{array}$$

leading to 16 periodic points in the two dimensional map, classified as six prime orbits

$$\begin{array}{c}00\\ 01\end{array},\begin{array}{c}00\\ 11\end{array},\begin{array}{c}10\\ 11\end{array},\begin{array}{c}01\\ 10\end{array},\begin{array}{c}00\\ 10\end{array},\begin{array}{c}01\\ 11\end{array}$$

their cyclic permutations

$$\begin{array}{c}01\\ 00\end{array},\begin{array}{c}11\\ 00\end{array},\begin{array}{c}11\\ 10\end{array},\begin{array}{c}10\\ 01\end{array},\begin{array}{c}10\\ 00\end{array},\begin{array}{c}11\\ 01\end{array}$$

and four repeated period one orbits

$$\begin{array}{c}00\\ 00\end{array},\begin{array}{c}01\\ 01\end{array},\begin{array}{c}10\\ 10\end{array},\begin{array}{c}11\\ 11\end{array}$$

### 2.3 Symmetry

Until now we have not considered the symmetry of the system. This is particularly important for cycle expansions , for which a correct treatment of symmetry significantly accelerates convergence of the expansions. An uncoupled lattice has symmetry group $`S_X\times G^X`$ where $`X`$ is the number of lattice sites (equal to two in this paper), $`S_X`$ is the symmetric group, of order $`X!`$, corresponding to permutations of the lattice points, and $`G`$ is the symmetry group of the dynamics. A symmetrically coupled lattice with $`X>2`$ has symmetry group $`D_X`$, the dihedral group, of order $`2X`$, corresponding to spatial cyclic permutations and reflections of the lattice sites. In the case $`X=2`$ the cyclic shift and reflection are equivalent, so the group is $`C_2`$, the cyclic group. Remnants of the internal symmetry group $`G`$ may or may not survive, depending on the nature of the coupling. For $`X>3`$ the coupling breaks the permutation symmetry, reducing the order of the symmetry group (ignoring internal symmetries), from $`X!`$ to $`2X`$ .

Each symmetry operation permutes periodic points. If the image of a periodic point lies in the same periodic orbit, that orbit is called symmetric, and is invariant under a cyclic subgroup of symmetry operations. If it does not, the new orbit is equivalent to the first under the symmetry, and shares many of the same properties, in particular stability and averages invariant under the symmetry. It is also possible for a periodic point to be a fixed point of the symmetry transformation, in which case its orbit is called a boundary orbit. In the presence of symmetry, the best strategy is desymmetrisation, in which all images of a phase space point are considered equivalent, and so may be described by a representative point on a subset of the original phase space called the fundamental domain. Equivalent orbits in the full dynamics correspond to a single orbit in the desymmetrised dynamics, and symmetric orbits in the full dynamics correspond to repeats of shorter orbits in the desymmetrised dynamics. Thus the enumeration of prime orbits up to a given period in the desymmetrised dynamics contains much more information than the enumeration up to the same period in the original dynamics, hence the accelerated convergence of cycle expansions.

We continue with the previous example, considering the $`N=2`$ maps, which have no internal symmetry. Thus the symmetry group is the cyclic group $`C_2`$, with the nontrivial element corresponding to $`(x,y)(y,x)`$. The fundamental domain can be taken as $`xy`$, from which the full domain is obtained by taking symmetry operations. The boundary satisfies $`x=y`$. Now the prime orbits of period two in the desymmetrised dynamics are classified as follows. There are two asymmetric orbits of period two in the original dynamics,

$$\begin{array}{c}00\\ 10\end{array},\begin{array}{c}10\\ 11\end{array}$$

one boundary orbit of original period two

$$\begin{array}{c}00\\ 11\end{array}$$

and three symmetric orbits of original period four

$$\begin{array}{c}00\\ 01\\ 00\\ 10\end{array},\begin{array}{c}01\\ 01\\ 10\\ 10\end{array},\begin{array}{c}10\\ 11\\ 10\\ 11\end{array}$$

Note that the total number (six) is the same, but the information content is higher, since all the original period two orbits are included as period one or two, and also some orbits of original period four.

In the case $`N=3`$, the symmetry group of the interacting case is $`C_2^2`$, since there is an internal symmetry $`xx`$, which in the coupled map survives as $`(x,y)(x,y)`$. This is in addition to the permutation $`(x,y)(y,x)`$. Thus the fundamental domain can be taken as $`x|y|`$. Orbits may now be classified as asymmetric, symmetric under $`(x,y)(y,x)`$, symmetric under $`(x,y)(x,y)`$, symmetric under $`(x,y)(y,x)`$, boundary with $`x=y`$, boundary with $`x=y`$ and (if it occurs) boundary with $`x=y=0`$.

### 2.4 Behaviour in the coupled case

In this section we illustrate the behaviour of the $`T_N`$ map periodic solutions when the coupling parameter $`a`$ is in the range $`[0,0.1]`$. In particular, we describe two phenomena, that is the way and stability and number of cycles change with $`a`$.

Stability is expressed by the real parts222henceforth we call $`\lambda `$ eigenvalue, but we always mean its real part. $`\lambda _1`$,$`\lambda _2`$ of the eigenvalues of the product of the Jacobian matrices

$$J(x)=\left(\begin{array}{cc}\frac{f(x,y)}{x}& \frac{f(x,y)}{y}\\ \frac{f(y,x)}{x}& \frac{f(y,x)}{y}\end{array}\right)$$

of each point of one cycle, where $`f(x,y)`$ and $`f(y,x)`$ are the set of equation (5,6). $`\lambda `$ is said to be stable when $`|\lambda |<1`$. When $`a=0`$, the systems is uncoupled and the $`T_N`$ maps have no stable cycles. Since the eigenvalues are continuous functions of $`a`$, we expect no stable eigenvalue to show up for very small $`a`$ at fixed orbit period. As for stability, we analyse the trend of the maximum and minimum eigenvalues found for each period. In general, we consider the product

$$\mathrm{\Lambda }=\underset{|\lambda _i|>1}{}\lambda _i$$

so that in the 2-d case it is either $`\lambda _1\lambda _2`$ or just the unstable $`\lambda `$. Stable orbits are not used for cycle expansions, hence they are not considered for this analysis. The graphs in Figure 2 show that both the maximum and the minimum eigenvalues increase hyperbolically (i.e. exponentially) with respect to the period $`n_p`$ of the cycles when $`a=0`$ and still up to $`a10^7`$, but the trend of the minimum $`\lambda `$’s changes as long as $`a`$ grows, so that their increase gets slower and slower. Besides, as we will explain later in this section, the few cycles with $`|\mathrm{\Lambda }|=4^{2n_p}`$ move off the fundamental domain when $`a0`$, so that the maximum eigenvalues take the form $`|\mathrm{\Lambda }|4^{n_p}`$.333The same happens with the $`T_3`$ coupled maps where the maximum $`|\mathrm{\Lambda }|`$ is $`9^{2n_p}`$ when $`a=0`$ and $`9^{n_p}`$ or lower otherwise. When $`a=0.1`$, the minimum eigenvalues are perhaps best fitted by a power law, rather than an exponential. Such a behaviour is common to all the coupled $`T_2`$ and $`T_3`$ maps. Figure 2 also shows that the trend of the minimum $`\lambda `$’s also gets less and less smooth, as $`a`$ increases. This phenomenon, together with the loss of hyperbolicity444for a hyperbolic system $`|\mathrm{\Lambda }|>Ae^{bn_p}`$ for some $`A,b>0`$, dramatically slows down the convergence of cycle expansions, as we discuss in the next section.

Figure 3 shows the trend of the eigenvalues in the $`2B^{}`$ map when $`a=0.1`$.

As we can see, it is very difficult to make a good fit for the trend of the minimum $`|\mathrm{\Lambda }|`$’s in this case, because all the coupled $`T_2`$ (and $`T_3`$) maps have got almost stable cycles for such $`a`$’s, including some at high period. It is interesting to compare this behaviour with what happens with intermittent $`1D`$ maps .

The other interesting issue we analyse is whether and why the number of cycles of a coupled map for a fixed period varies as a function of the coupling parameter $`a`$. Numerical simulation shows that few periodic orbits of the uncoupled map (7) move off the fundamental domain when the system is coupled ($`a>0`$). Moreover, the number of $`n_p4`$ cycles continues to decrease as $`a`$ grows, because some periodic solutions of the set (5,6) become complex or degenerate for certain values of $`a`$, so that most of the missing orbits have just disappeared into the complex plane. Data representing the number of periodic points $`H(a,n_p)`$ with $`a[10^6,10^3]`$ have been fitted with the curve

$$H(a,n_p)=N^{2n_p}\mathrm{exp}(2n_pCa^\eta )$$
(8)
where $`C`$ and $`\eta `$ are parameters that depend on the model employed, as shown in table 1.

As we can see, $`\eta =0.5`$ in all the models, which is theoretically plausible since the dynamics and bifurcations are controlled by the critical point of the Tchebyscheff map, which is quadratic. That implies that perturbations of order $`a`$ have effects of order $`\sqrt{a}`$, hence the exponent of $`0.5`$. $`H(a,n_p)`$ is related to the topological entropy $`h(a)`$ of the system by

$$h(a)=\underset{n_p\mathrm{}}{lim}\frac{\mathrm{ln}[H(a,n_p)]}{n_p}$$

so that

$$h(a)=2[\mathrm{ln}NCa^\eta ]$$

which decreases with a power law as $`a`$ increases. Figure 4 shows good agreement of the $`h(a)`$ fitting function with the data for very weak coupling ($`a10^3`$) while a slight deviation for larger $`a`$. This deviation is not unexpected, since this form for the topological entropy is only an approximate expression in the small coupling limit. Of course a detailed investigation could reveal that the topological entropy is not smooth in fine detail, however this form is remarkably accurate for practical purposes.

## 3 cycle expansions and computation of averages

### 3.1 cycle expansions

Given a quantity $`v(x)`$, we want to compute its average on the attractor of the dynamics (5,6), such as the one depicted in figure 5. Assuming that such attractor is ergodic for any value of the coupling parameter we consider, direct simulation is the quickest method to use for this purpose, but its precision is limited by statistics. That is one reason why we use cycle expansions, which are potentially far more precise. According to periodic orbit theory , the average of $`v`$ can be found by computing the dynamical zeta function

$$1/\zeta =\underset{p}{}\left(1t_p\right),t_p=\frac{z^{n_p}e^{\beta V_p}}{|\mathrm{\Lambda }_p|}$$
(9)
where

$$V_p=\underset{i=0}{\overset{n1}{}}v(f^i(x_0))$$
(10)
$`f`$ is the map in question, $`\beta `$ is a parameter, $`\mathrm{\Lambda }_p`$ the eigenvalue of the cycle $`p`$ and $`n`$ is the orbit period. The uncertain digit in the average is obtained by comparing results from the computation of (9) up to the two largest periods (for example 10 and 11). It is possible to check the rate of convergence of the zeta function by setting $`\beta =0`$ and $`z=1`$ and computing $`1/\zeta (0,0)`$, which should approach zero. This way we can have an idea of how effective our cycle expansion is and whether longer cycles are needed. Figure 6 shows the convergence of $`1/\zeta (0,0)`$ versus the truncation period of the expansion.

As we can see, the convergence is monotonic and exponential when the coupling parameter is relatively small (say $`a10^4`$), while it becomes nonmonotonic and increasingly slower for larger values of $`a`$.

The accuracy of $`1/\zeta `$ is particularly affected towards $`a=0.1`$, where longer cycles or stability ordering may be needed. Ordering by stability means to expand the zeta function keeping all terms whose period and stability are less than a given cutoff ($`|\mathrm{\Lambda }|<\mathrm{\Lambda }_{max}`$). This scheme is usually preferred to length ordering in case of intermittency or bad symbolic dynamics grammar , but it hasn’t been attempted here since it would need much a higher period of truncation for the cycle expansion and, as a consequence, impractical CPU times. A full resolution would require understanding which orbits are close to stability, which is beyond the scope of this paper.

Anyway the behaviour shown in figure 6 traces back to the trend of the minimum eigenvalues of the cycles (sec. 2.3), since we can see from equation (9) that the smaller $`|\mathrm{\Lambda }|`$’s are more relevant in the computation of the zeta function. That means that if their trend is exponential with respect to increasing $`a`$’s, the $`t_p`$’s are smaller and smaller as $`1/\zeta (0,0)`$ approaches zero. Otherwise the convergence gets more complicated. The expansion for the $`T_3`$ maps has been truncated at period $`7`$, yet the precision obtained is comparable with what we have got in the $`T_2`$ maps. That is because $`9^74^{11}`$, meaning a similar number of periodic points and hence more or less the same number of terms in the two expansions.

### 3.2 Computation of averages

We now want to use cycle expansions to evaluate the average of an observable over the two-site lattice (5,6) and compare the outcomes with what we get with direct simulation on the same lattice. According to the theory of Beck , an interaction potential

$$W(a)=\frac{1}{2}\mathrm{\Phi }^i\mathrm{\Phi }^{i+1}$$
(11)
and self-interacting potentials

$`V^{(N=2)}(a)={\displaystyle \frac{2}{3}}\mathrm{\Phi }^3+\mathrm{\Phi }`$ (12)
$`V^{(N=3)}(a)=\mathrm{\Phi }^4+{\displaystyle \frac{3}{2}}\mathrm{\Phi }^2`$ (13)
can be defined for the chaotic strings described by equation (4). Zeros of (11) and minima of (12,13) may represent parameters of the Standard Model of particle physics such as coupling constants and masses of elementary particles. Since both the potentials defined above depend on $`a`$ in a nontrivial way, we can only find zeros and minima numerically by averaging both $`V(a)`$ and $`W(a)`$ over the lattice (4). Beck has obtained his results by running a long trajectory on a lattice with many sites, we are going to see what happens instead by computing the zeta function (9) on the two-site lattice (5,6).

Figure 7 shows our results (achieved by computing (9) up to two different orbit periods) compared with direct simulation of the potentials $`V^{(N)}`$.A similar behaviour has been observed with the interaction potential $`W(a)`$. Five main features emerge.

1. In most points cycle expansion coincides with direct simulation.
2. Cycle expansion is far offset on some points, which can be explained with the presence of almost stable cycles ($`|\mathrm{\Lambda }|1`$) of high periods whose terms in eq. (9) have too much weight in the expansion.
3. There are points where both cycle expansion and direct simulation are offset. That is probably caused by one or more stable cycles or chaotic attractors, whose basins of attraction alter the phase space making it no longer ergodic.
4. The precision of the expansion visibly decreases as $`a`$ approaches $`0.1`$, as we expect from the observations made in the previous sections on the trend of the minimum eigenvalues and the behaviour of $`1/\zeta (0,0)`$.
5. Figure 8 shows an offset at some point between cycle expansion and direct simulation in the model $`2B^+`$ performed with different starting values from those used in figure 7. That proves a lack of ergodicity in some regions of the interval $`[0,0.1]`$ of $`a`$ and indicates that the unstable periodic orbits are mostly unaffected by the presence of a stable orbit (or other ergodic component such as an elliptic region). Hence the cycles form a transient chaotic set more or less covering the attractor (Fig. 5) for similar values of the parameter. In this sense, the periodic orbit expansion results are more stable than direct simulation with respect to variation of the parameter, and may be preferred for this reason.

## 4 Conclusions and outlook

Periodic orbits have provided a unique window with which to view properties of weakly coupled CML. The periodic orbits in the uncoupled limit provide a reliable means of finding periodic orbits of the weakly coupled system. The symmetry of the coupled system has been exploited to improve convergence of the cycle expansion. Cycle expansion techniques work as expected, well near the uncoupled, more chaotic limit, but with slower convergence at stronger coupling. The computation of averages using cycle expansions leads to important information at a glance concerning the number and nature of nearly or fully stable states coexisting with the chaotic web of unstable periodic orbits.

The observed power law dependence of the topological entropy with parameter requires a more detailed investigation of the density and nature of bifurcations in this system to explain the value of the power $`\eta =1/2`$, the extent of its apparent universality and corrections for finite values of the parameter. For a more rigorous approach to the topological entropy of a different two site model, see .

This work has clarified a number of issues and challenges regarding a more general periodic orbit approach to CML. For stronger coupling (hence weaker chaos) stability ordering of cycle expansions  is likely to be more effective than the length ordering employed here; this requires an understanding of the symbolic dynamics of the most stable cycles, which may be obtained in the study of the bifurcations mentioned above. Larger lattices are of course an important aspect, however the rapid decay of spatial correlations for weakly coupled CML make this less important than it might seem. The fact, as observed here, that small windows of stability at weak coupling affect the cycle expansions at only a few parameter values, gives great confidence since their presence is expected to be even less evident in lattices of larger size, in which typical orbits have many unstable directions.

## 5 Acknowledgements

We acknowledge financial support of the Nuffield Foundation (grant NAL/00353/G) and helpful discussions with P. Cvitanović. DL would also like to thank G. Servizi and V. Smith for having made available computing facilities in Bologna and Bristol respectively, as well as N. Chevalier, R. Morris and A. Shall for technical support.


text:
# 𝐵→ϕ𝐾_𝑆 𝐶𝑃 asymmetries as a probe of supersymmetry1footnote 11footnote 1Talk presented at ICFP2003, Seoul, Oct 6–11, 2003, to appear in the proceedings.

## Abstract

We consider sparticle contributions to time dependent $`CP`$ asymmetry in $`B\varphi K_S`$. As for the gluino-squark loop, $`LR`$ or $`RL`$ insertion is more likely to give recently observed $`S_{\varphi K}<0`$ than $`LL`$ or $`RR`$ insertion. Neutral Higgs contribution does not change $`S_{\varphi K}`$ very much due to the strong constraint from $`B_s\mu ^+\mu ^{}`$. We also show correlations among related processes such as direct $`CP`$ asymmetry in $`BX_s\gamma `$ and $`B_s`$$`\overline{B_s}`$ mixing, and discuss theoretical motivations for desired strength of flavor violation.

preprint: KAIST-TH 2003/12

Time dependent $`CP`$ asymmetry in $`B\varphi K_S`$ mode is written as

$`𝒜_{\varphi K}(t)`$ $`{\displaystyle \frac{\mathrm{\Gamma }(\overline{B}^0(t)\varphi K_S)\mathrm{\Gamma }(B^0(t)\varphi K_S)}{\mathrm{\Gamma }(\overline{B}^0(t)\varphi K_S)+\mathrm{\Gamma }(B^0(t)\varphi K_S)}}`$
$`=C_{\varphi K}\mathrm{cos}(\mathrm{\Delta }m_dt)+S_{\varphi K}\mathrm{sin}(\mathrm{\Delta }m_dt),`$
where $`\mathrm{\Delta }m_d`$ denotes the mass splitting between the two mass eigenstates composed of $`B^0`$ and $`\overline{B}^0`$. Current data on the coefficients of the sine and cosine terms and their Standard Model (SM) predictions are summarized in Table 1.

It shows a 2.7 $`\sigma `$ discrepancy between the measured value and SM prediction of $`S_{\varphi K}`$. This has become a hot issue in the high energy physics community recently. One of the reasons why physicists are focusing on $`B\varphi K_S`$ mode is that it is loop suppressed in the SM because it is a purely flavor changing neutral current (FCNC) process. This means that it is more sensitive to possible new physics effects which presumably arise at loop level, than those processes that have tree level contribution such as $`BJ/\psi K_S`$. In this work, we consider a possibility of filling this gap between the data and SM value, with effects from supersymmetry (SUSY).

One major source of flavor violation in the general minimal supersymmetric standard model (MSSM) is non-diagonal squark mass matrix. An off-diagonal element in the squark mass matrix in the super CKM basis can constitute a diagram leading to $`bs\overline{q}q`$ through gluino mediation, where $`q=u,d,s,c,b`$. An example of this type of QCD penguin diagram is shown in Fig. 1.

Here we adopt the mass insertion (MI) approximation which is one of the most convenient ways to analyzing this type of contribution. We consider four MI’s relevant to $`bs`$ transitions, i.e., $`(\delta _{LL}^d)_{23}`$, $`(\delta _{RR}^d)_{23}`$, $`(\delta _{LR}^d)_{23}`$, and $`(\delta _{RL}^d)_{23}`$, one at a time. In addition to these four cases, we include another one called $`RL`$ dominance scenario Everett:2001yy . In this case, $`C_{7\gamma }`$ and $`C_{8g}`$ are assumed to vanish at the $`m_b`$ scale, so the observed decay of $`BX_s\gamma `$ must come from $`\stackrel{~}{C}_{7\gamma }`$, the chirality-flipped version of $`C_{7\gamma }`$. This unconventional scenario and the usual case can be distinguished by measuring the photon polarization. Therefore $`(\delta _{RL}^d)_{23}`$ is forced to be a finite size to account for the data unlike the other four cases, and this corresponds to an extreme case with maximized SUSY effect in some sense. In the following, we turn on one of the four MI parameters at a time, and scan over it imposing constraints from $`BX_s\gamma `$ and $`B_s`$$`\overline{B_s}`$ mixing as follows:

$`2.0\times 10^4<B(BX_s\gamma )<`$ $`4.5\times 10^4,`$
$`\mathrm{\Delta }M_s>14.9\text{ps}^1`$ $`\text{Stocchi:2000ps}.`$
We impose a rather generous bound on $`B(BX_s\gamma )`$ to take into account theoretical uncertainties. The common squark mass $`\stackrel{~}{m}`$ and the gluino mass $`m_{\stackrel{~}{g}}`$ are chosen to be $`m_{\stackrel{~}{g}}=\stackrel{~}{m}=400`$ GeV. We use QCD factorization Beneke:1999br  in evaluating hadronic matrix elements for $`B\varphi K_S`$. We do not assume any new physics in the $`B^0`$$`\overline{B}^0`$ mixing because $`\mathrm{sin}2\beta `$ measurement from the $`BJ/\psi K_S`$ mode is well consistent with the SM fit.

First, let us consider the $`LL`$ insertion case. Fig. 2 (a) is the region of $`(\delta _{LL}^d)_{23}`$ consistent with $`B(BX_s\gamma )`$ and $`\mathrm{\Delta }M_s`$ constraints. The solid lines are contours of $`S_{\varphi K}`$. The small square at the origin corresponds to the SM case, because there is no SUSY contribution there. The white space on the left is excluded by too small value of $`B(BX_s\gamma )`$. The two holes around the $`\mathrm{Re}(\delta _{LL}^d)_{23}=0`$ line are excluded by too small $`\mathrm{\Delta }M_s`$. In Fig. 2 (b), we plot predictions of $`S_{\varphi K}`$ and $`C_{\varphi K}`$ from the allowed values of $`(\delta _{LL}^d)_{23}`$. The dotted horizontal line is the current 1-$`\sigma `$ upper limit on $`S_{\varphi K}`$. This figure shows that the minimal value of $`S_{\varphi K}`$ in this case is about +0.5, and the change of $`S_{\varphi K}`$ from the SM prediction is not significant. Although $`S_{\varphi K}`$ is not very much affected, large effect in $`B_s`$$`\overline{B_s}`$ can be expected. For example, we show the prediction of $`\mathrm{\Delta }M_s`$ in Fig. 2 (c). The SM value of $`\mathrm{\Delta }M_s`$ is about 16 $`\mathrm{ps}^1`$, marked by the black square. We can read huge enhancement of $`\mathrm{\Delta }M_s50\mathrm{ps}^1`$ is possible. Also, we show $`\mathrm{sin}2\beta _s`$ in Fig. 2 (d), where $`2\beta _s`$ is the phase of $`B_s|H_{\mathrm{eff}}^{\mathrm{\Delta }B=2}|\overline{B_s}`$. The SM prediction of it is almost zero, but it can have any value between $`1`$ and 1 here, which means that $`CP`$ violation in $`B_s`$$`\overline{B_s}`$ mixing can drastically change.

The $`RR`$ insertion case is almost identical to the $`LL`$ insertion case. The only difference is the way $`B(BX_s\gamma )`$ constrains $`(\delta _{RR}^d)_{23}`$. The $`LL(RR)`$ insertion mainly contributes to $`C_{7\gamma }(\stackrel{~}{C}_{7\gamma })`$, and the dependence is such that

$`B(BX_s\gamma )`$ $`|C_{7\gamma }^{\mathrm{SM}}+C_{7\gamma }^{\mathrm{SUSY}}|^2,\text{for }LL,`$
$`B(BX_s\gamma )`$ $`|C_{7\gamma }^{\text{SM}}|^2+|\stackrel{~}{C}_{7\gamma }^{\text{SUSY}}|^2,\text{for }RR.`$
Because of this difference, $`(\delta _{RR}^d)_{23}`$ plane does not show a region excluded by $`B(BX_s\gamma )`$ that is present in Fig. 2 (a). However, this additional parameter space does not help shifting down $`S_{\varphi K}`$ very much. Plots for this case are available in Ref. Kane:2002sp .

It is well known that the chirality-flipping $`LR`$ and $`RL`$ insertions receive enhancement by the factor of $`m_{\stackrel{~}{g}}/m_b`$, when they contribute to a (chromo) magnetic dipole operator. This is why they are strongly constrained by $`BX_s\gamma `$. For the same reason, they are more effective in modifying $`S_{\varphi K}`$ and/or $`C_{\varphi K}`$ than $`LL`$ or $`RR`$ insertion given the same magnitude of the MI parameter. We show the region of $`(\delta _{LR}^d)_{23}`$ consistent with $`B(BX_s\gamma )`$ and $`\mathrm{\Delta }M_s`$ in Fig. 3 (a). As in the previous case, the solid lines are contours of $`S_{\varphi K}`$, and the small square at the origin corresponds to the SM case. The annulus with radius $`10^2`$ comes from the $`B(BX_s\gamma )`$. Here the $`B(BX_s\gamma )`$ constraint is so strong that $`\mathrm{\Delta }M_s`$ does not change very much from the SM value and play no role in constraining $`(\delta _{LR}^d)_{23}`$. We show $`S_{\varphi K}`$ and $`C_{\varphi K}`$ from these values of $`(\delta _{LR}^d)_{23}`$ in Fig. 3 (b). The dotted box is 1-$`\sigma `$ intervals of these observables, and we find some region where both $`S_{\varphi K}`$ and $`C_{\varphi K}`$ are within the box. This plot also shows a definite correlation between them. For $`S_{\varphi K}`$ less (bigger) than the SM value, $`C_{\varphi K}<(>)0`$. Fig. 3 (c) is the plot of $`B(B\varphi K_S)`$, with hatches on the region where $`B(B\varphi K_S)>1.6\times 10^5`$. This region, excluded by $`B(B\varphi K_S)`$ constraint, corresponds to the hatched region in Fig. 3 (a). Here we learn that significant portion of the $`(\delta _{LR}^d)_{23}`$ parameter space that is consistent with $`B(BX_s\gamma )`$, is excluded by $`B(B\varphi K_S)`$. In Fig. 3 (d) is displayed the correlation between $`S_{\varphi K}`$ and direct $`CP`$ asymmetry in $`BX_s\gamma `$. For $`S_{\varphi K}`$ less (bigger) than the SM value, $`A_{\mathrm{CP}}^{bs\gamma }>(<)0`$. It appears that $`A_{\mathrm{CP}}^{bs\gamma }`$ is outside the dotted box, but the experimental uncertainty is still large, so we cannot definitely conclude that this scenario is disfavored currently.

We turn to the $`RL`$ dominance scenario. Fig. 4 (a) is the region of $`(\delta _{RL}^d)_{23}`$ constrained by $`B(BX_s\gamma )`$ and $`\mathrm{\Delta }M_s`$. Again, the solid lines are contours of $`S_{\varphi K}`$. In this case, we are assuming that SM contributions to $`C_{7\gamma }`$ and $`C_{8g}`$ are somehow canceled by SUSY contributions. Because we are always involving SUSY contributions in $`C_{7\gamma }`$ and $`C_{8g}`$, there cannot be any point on the $`(\delta _{RL}^d)_{23}`$ plane that reduces to the SM case. The allowed annulus, centered at the origin, has radius $`10^2`$, which is fixed by $`B(BX_s\gamma )`$. As in the previous case, the region consistent with $`B(BX_s\gamma )`$ is always consistent with $`\mathrm{\Delta }M_s`$ as well. Predictions of $`S_{\varphi K}`$ and $`C_{\varphi K}`$ from the allowed region are shown in Fig. 4 (b). We can fit $`S_{\varphi K}`$ to any value between $`1`$ and 1. Moreover, big change in $`C_{\varphi K}`$ from the SM value of zero, is expected in general, resulting in any value between $`1`$ and 1. However, some region around the center is not covered. In Fig. 4 (c), we plot $`B(B\varphi K_S)`$, with hatches on the region with excessive $`B(B\varphi K_S)`$. Since there is only one weak phase in the $`bs\gamma `$ amplitude, we have $`A_{\mathrm{CP}}^{bs\gamma }=0`$ in this case. If we turned on another weak phase such as coming from $`(\delta _{RR}^d)_{23}`$, we could have non-vanishing $`A_{\mathrm{CP}}^{bs\gamma }`$, but we do not pursue more detailed analysis here.

We also analyzed the usual $`RL`$ insertion case, keeping SM contributions to $`C_{7\gamma }`$ and $`C_{8g}`$. Because the SM prediction of $`B(BX_s\gamma )`$ is already consistent with the data, there is less room for new physics in this case than in the $`RL`$ dominance scenario. Nevertheless, we can get sizable change in $`S_{\varphi K}`$, satisfying $`B(BX_s\gamma )`$. Plots for this case can be found in Ref. Kane:2002sp .

These results of analysis depend on the sparticle mass scale. For a fixed value of MI parameter, SUSY effects get enhanced as the sparticle masses decrease. In Figs. 5, we show possible values of $`S_{\varphi K}`$ as functions of the gluino mass, $`m_{\stackrel{~}{g}}`$.

In these plots, we fix the ratio $`xm_{\stackrel{~}{g}}^2/\stackrel{~}{m}^2`$ to be 1, and scan over one of the four MI parameters imposing constraints from $`B(BX_s\gamma )`$ and $`\mathrm{\Delta }M_s`$. Figs. 5 (a) and (b) are for the chirality-preserving insertions, $`LL`$ and $`RR`$, respectively. Here we see a clear tendency of SUSY effects decoupling as the gluino mass increases. Note that our analysis shown in Figs. 24 is for $`m_{\stackrel{~}{g}}=400`$ GeV. The way to getting $`S_{\varphi K}`$ consistent with the data is to imagine the gluino mass $`250`$ GeV, which is close to the current experimental lower bound. On the other hand, a chirality-flipping insertion, $`LR`$ or $`RL`$, shown in Figs. 5 (c) and (d), results in a range of $`S_{\varphi K}`$ that is independent of $`m_{\stackrel{~}{g}}`$. This is explained in the following way in the $`LR`$ insertion case for example. The dependence of $`BX_s\gamma `$ on the SUSY parameters is only through the ratio of $`(\delta _{LR}^d)_{23}/m_{\stackrel{~}{g}}`$, which precisely is the parameter that $`S_{\varphi K}`$ also depends on. The set of $`(\delta _{LR}^d)_{23}/m_{\stackrel{~}{g}}`$ allowed by $`B(BX_s\gamma )`$ is independent of $`m_{\stackrel{~}{g}}`$, and so is $`S_{\varphi K}`$. The $`RL`$ and $`RL`$ dominance cases are explained likewise. In fact, this type of constant dependence would be true for $`LL`$ or $`RR`$ insertion case as well, if we allowed arbitrarily large size of $`(\delta _{LL}^d)_{23}`$ or $`(\delta _{RR}^d)_{23}`$ and if we discarded $`\mathrm{\Delta }M_s`$ constraint.

In the QCD factorization, divergences coming from endpoint of integral $`𝑑x/x`$ are regulated in the form of $`X_{A(H)}=(1+\varrho _{A(H)}e^{i\phi _{A(H)}})\mathrm{ln}(m_B/\mathrm{\Lambda }_h)`$, where $`\mathrm{\Lambda }_h0.5`$ GeV is the scale in hard-scattering. Discussion so far has been carried out without considering any uncertainties coming from this regularization, but we have to scan over the parameters $`\varrho _{A,H}`$ and $`\phi _{A,H}`$ properly to take care of these theoretical uncertainties. In Fig. 6 (a), we show a curve of $`S_{\varphi K}`$ as a function of the gluino mass. In this plot, we fix $`x=1`$, and turn on only the $`RR`$ insertion with the value $`(\delta _{23}^d)_{RR}=0.5340.856i`$, which minimizes $`S_{\varphi K}`$ at $`m_{\stackrel{~}{g}}=200`$ GeV. This curve is broadened into bands shown in Fig. 6 (b), if we begin to scan over the aforementioned parameters. In general, the uncertainty in the annihilation from $`\varrho _A`$ dominates over that in hard spectator scattering from $`\varrho _H`$. The black band is for the range $`0\varrho _A,\varrho _H1`$, and this is what the authors of Ref. Beneke:1999br  recommend to use. On this band, the tendency of decoupling SUSY effects with increasing $`m_{\stackrel{~}{g}}`$ is still maintained. However, for the case of $`0\varrho _A,\varrho _H5`$, which is displayed in light gray, estimated theoretical uncertainties are so large that almost every value of $`S_{\varphi K}`$ between $`1`$ and 1 is predicted, for any gluino mass between 200 GeV and 600 GeV. Ref. Ciuchini:2002uv  claims that they get $`S_{\varphi K}<0`$ from $`LL`$ or $`RR`$ insertion with $`m_{\stackrel{~}{g}}=350`$ GeV, disagreeing with us. We suspect this is due to the difference in choosing intervals for $`\varrho _A`$ and $`\varrho _H`$. They use $`0\varrho _A8`$, and from Fig. 6 (b), it is obvious that $`S_{\varphi K}`$ will have every value between $`1`$ and 1 for this interval.

Let us think of theoretical motivations for $`(\delta _{LR}^d)_{23}`$ or $`(\delta _{LR}^d)_{23}10^2`$. In many alignment models using flavor symmetry, there remain small off-diagonal elements in the squark mass matrix in the super CKM basis. They can cause $`(\delta _{23}^d)_{LL}`$ or $`(\delta _{23}^d)_{RR}`$ $`10^2`$ with arbitrary complex phase. This chirality preserving insertion can lead to an induced $`LR`$ or $`RL`$ insertion,

$$(\delta _{LR}^d)_{23}^{\mathrm{ind}}=(\delta _{LL}^d)_{23}\times \frac{m_b(A_b\mu \mathrm{tan}\beta )}{\stackrel{~}{m}^2}10^2,$$

provided $`\mathrm{tan}\beta `$ is large enough so that $`\mu \mathrm{tan}\beta 30`$ TeV. This kind of double insertion mechanism has been used in explaining both $`ϵ_K`$ and $`ϵ^{}/ϵ_K`$ from a single $`CP`$ violation source of $`(\delta _{LL}^d)_{12}`$ Baek:1999jq . We also constructed a model that naturally gives the desired amount of $`LR`$ or $`RL`$ insertion using intersecting D5 branes Kane:2002sp .

Since QCD penguin contributes to $`bs\overline{q}q`$ for $`q=u,d,s,c,b`$ at equal strength, one may worry about what happens to other decays such as $`B\eta ^{}K_S,\pi K,K^+K^{}K_S`$. There are other 4-quark operators which contribute only to these decays not affecting $`B\varphi K_S`$, and changes in these decay modes are not definitely predictable in the present study. On the other hand, there is a nice mechanism which enables us to control $`BVP`$ and $`BPP`$ modes independently. Parity invariance tells us that the transition amplitudes of these modes depend on the new physics contribution in such a way that $`BVP(PP)`$ mode depends on $`_i[C_i^{\mathrm{NP}}+()\stackrel{~}{C}_i^{\mathrm{NP}}]\mathrm{\#}_i`$, where $`C_i^{\mathrm{NP}}`$ and $`\stackrel{~}{C}_i^{\mathrm{NP}}`$ are Wilson coefficients with opposite chiralities, coming from new physics. Hence we can control $`B\eta ^{}K_S`$ and $`B\varphi K_S`$ separately, and only the latter will change if we imagine a situation where $`(\delta _{LR}^d)_{23}=(\delta _{RL}^d)_{23}`$, for instance. If this is the case, predictions of other observables such as $`C_{\varphi K}`$ and $`A_{\mathrm{CP}}^{bs\gamma }`$ will be different than the four cases we considered here. After all, it requires further study to see changes in other decay modes in more general cases.

From now, we consider the neutral Higgs mediated contribution to $`B\varphi K_S`$. A typical diagram is shown in Fig. 7.

This diagram contributes to $`B_s\mu ^+\mu ^{}`$ as well as to $`B\varphi K_S`$ if we replace the $`s\overline{s}`$ pair with $`\mu ^+\mu ^{}`$. Once we impose the upper bound from CDF Abe:1998ah , $`B(B_s\mu ^+\mu ^{})<2.6\times 10^6`$, we find that $`S_{\varphi K}>0.71`$. This argument applies not only in minimal flavor violation scenario, but also in general MSSM Kane:2002sp . As a result, Higgs exchange does not cause substantial change in $`S_{\varphi K}`$.

In conclusion, we analyzed SUSY contributions to time dependent $`CP`$ asymmetry in $`B\varphi K_S`$. As to the gluino-squark loop, negative $`S_{\varphi K}`$ is more likely to come from $`LR`$ or $`RL`$ insertion than from $`LL`$ or $`RR`$ insertion if the sparticle masses are not close to the current experimental lower bound. However, we may have chances to look for $`LL`$ or $`RR`$ insertion in $`B_s`$$`\overline{B_s}`$ mixing. Correlations among $`S_{\varphi K}`$, $`C_{\varphi K}`$, $`A_{\mathrm{CP}}^{bs\gamma }`$, $`\mathrm{\Delta }M_s`$, $`\mathrm{sin}2\beta _s`$ may help us discriminate among different possibilities. Constraint from nonleptonic decays such as $`B\varphi K_S`$ is becoming as severe as that from $`BX_s\gamma `$. Higgs mediated FCNC cannot explain $`S_{\varphi K}<0`$. Finally, let us direct the reader to Ref. Kane:2002sp  for more complete discussion including dilepton $`CP`$ asymmetry in $`\mathrm{{\rm Y}}(4S)`$ decay.

The author is grateful to G. L. Kane, P. Ko, C. Kolda, Haibin Wang, and Lian-Tao Wang, for enjoyable collaboration. This work was supported in part by KOSEF through CHEP at Kyungpook National University and by the BK21 program.


text:
# Quantum key distribution using gaussian-modulated coherent states

#### Acknowledgements.

The contributions of J. Gao to early stages of the experiment, and of K. Nguyen to the software development, are acknowledged. We thank S. Iblisdir for discussions, and Th. Debuisschert for the loan of the 780 nm integrated modulator. This work was supported by the EU programme IST/FET/QIPC (projects “QUICOV” and “EQUIP”), the French programmes ACI Photonique and ASTRE, and by the Belgian programme ARC.

#### Correspondence

and requests for materials should be addressed to Philippe Grangier (e-mail philippe.grangier@iota.u-psud.fr).

## Appendix A Methods

### Relevant Heisenberg relations

In a RR protocol, Alice’s estimator for $`x_B`$ and Eve’s estimator for $`p_B`$ can be denoted respectively as $`\alpha x_A`$ and $`\beta p_E`$, where $`\alpha `$, $`\beta `$ are real numbers. The corresponding errors are $`x_{B|A,\alpha }=x_B\alpha x_A`$, and $`p_{B|E,\beta }=p_B\beta p_E`$. Because Alice’s, Bob’s, and Eve’s operators commute, we have $`[x_{B|A,\alpha },p_{B|E,\beta }]=[x_B,p_B]`$, and thus the Heisenberg relation $`\mathrm{\Delta }x_{B|A,\alpha }^2\mathrm{\Delta }p_{B|E,\beta }^2N_0^2`$. Defining the conditional variances as $`V(x_B|x_A)=\mathrm{min}_\alpha \{\mathrm{\Delta }x_{B|A,\alpha }^2\}`$ and $`V(p_B|p_E)=\mathrm{min}_\beta \{\mathrm{\Delta }p_{B|E,\beta }^2\}`$, we obtain $`V(x_B|x_A)V(p_B|p_E)N_0^2`$, or, by exchanging $`x`$ and $`p`$, $`V(p_B|p_A)V(x_B|x_E)N_0^2`$.

Alice has the estimators $`(x_A,p_A)`$ for the field $`(x_{\mathrm{in}},p_{\mathrm{in}})=(x_A+A_x,p_A+A_p)`$ that she sends, with $`A_x^2=A_p^2=sN_0`$. Here $`s`$ measures the amount of squeezing possibly used by Alice in her state preparation 14 , with $`sV^1`$ for consistency with Heisenberg’s relations. By calculating $`p_A^2=(Vs)N_0`$, $`p_B^2=G_p(V+\chi _p)N_0`$, $`p_Ap_B=G_p^{1/2}p_A^2`$, we obtain the conditional variance $`V(p_B|p_A)=p_B^2|p_Ap_B|^2/p_A^2=G_p(s+\chi _p)N_0`$. This equation and the constraint $`sV^1`$ gives $`V(p_B|p_A)G_p(V^1+\chi _p)N_0`$, and similarly $`V(x_B|x_A)G_x(V^1+\chi _x)N_0`$. The bound on $`V_{B|A}`$ is thus obtained by assuming that Alice may use squeezed or entangled beams, while the bound on $`V_{B|E}`$ can only be achieved if Eve uses an entangling attack. This reflects the fact that squeezing or entanglement play a crucial role in our security demonstration, even though the protocol implies coherent states. Our security proof addresses individual gaussian attacks only, but as the entangling cloner attack saturates the Heisenberg uncertainty relations, we conjecture that it encompasses all incoherent (non-collective) eavesdropping strategies.

### Experimental set-up

A continuous-wave laser diode at 780 nm wavelength associated with an acousto-optic modulator is used to emit 120-ns (full-width at half-maximum) pulses at a 800 kHz rate. The signal pulses contain up to 250 photons, while the local oscillator (LO) power is $`1.3\times 10^8`$ photons per pulse. The amplitude of each pulse is arbitrarily modulated by an integrated electro-optic modulator. However, owing to the unavailability of a fast phase modulator at 780 nm, the phase is not randomly modulated but scanned continuously. No genuine secret key can be distributed, strictly speaking, but random permutations of Bob’s data are used to provide realistic data (see Appendix B). The data are organized in bursts of 60,000 pulses, separated by synchronization periods also used to lock the phase of the LO. The overall homodyne detection efficiency is 0.81, due to the optical transmission (0.92), the mode-matching efficiency (0.96) and the photodiode quantum efficiency (0.92). For the critical data at 3 dB loss, the mode-matching efficiency was improved to 0.99, and thus the overall efficiency was 0.84. We also point out that many blocks of data were exchanged around the 3 dB loss point, with a typical rate above 55 kbit s-1.

### Secret key distillation

A common bit string is extracted from the continuous data by sequentially reconciling several strings (“slices”) of binary functions of the gaussian key elements, applying a binary reconciliation protocol successively on each bit 8 ; 10 . Here, we used five slices, each being corrected either by a trivial one-way protocol (communicating the bits) when the bit error rate (BER) is high, or by the two-way protocol Cascade 27 ; 28  when the BER is low. Note that the disclosed slices are useful for reconciling the remaining slices with less information leaking to Eve, even though they of course do not yield secret bits as such. In addition, Alice and Bob encrypt their classical messages using the one-time pad scheme with a fraction of the previous key bits, or a bootstrap key for the first block. For slices corrected with Cascade, the exchanged parities are encrypted with the same key bits on both sides 29 , making Eve aware of the differences between Alice’s and Bob’s parities (that is, the error positions) but not of their individual values. Fully communicated slices are also encrypted, thereby revealing no information at all to Eve. Still, Eve may exploit the interactivity of Cascade and gain some information on the final key by combining her knowledge of the error positions with that of the correlations between Alice’s and Bob’s gaussian values. In the present protocol, this information is numerically calculated for an entangling cloner attack, and then destroyed by privacy amplification. This is achieved by appropriate “hashing” 30  functions (see Appendix B). The resulting net secret key rate is then obtained by subtracting, from the raw key rate, the cost of the one-time pad encryption and the error-position information. Finally, let us emphasize that sliced reconciliation can be made very close to a one-way protocol by increasing the number of key elements from which the bits are jointly extracted (multidimensional reconciliation 8 ). This approach was not implemented here, but should deliver an improved secret key rate, approaching the value from the Csiszár-Körner formula 21 ; 22 .

## Appendix B Supplementary Information

### Experimental set-up

The source consists of a CW laser diode (SDL 5412) at 780 nm associated with an acousto-optic modulator, used to chop pulses with a duration 120 ns (full width at half-maximum), at a repetition rate 800 kHz. To reduce the excess noise, a grating-extended external cavity is used, and the beam is spatially filtered using a single mode fiber. Light pulses are then split onto a beam-splitter, one beam being the local oscillator (LO), the other Alice’s signal beam. The data is organised in bursts of 60,000 pulses, separated by time periods used to lock the phase of the LO, and sequences of pulses to synchronize the parties. In the present experiment, there is a burst every 1.6 seconds, corresponding to a duty cycle of about 5%, which is obviously under-optimised but should be easy to improve in further experiments.

The coherent state distribution is generated by modulating both the amplitude and phase of the light pulses with the appropriate probability law. In the present experiment, the amplitude of each pulse is arbitrarily modulated at the nominal 800 kHz rate by an integrated electro-optic LiNbO3 Mach-Zehnder interferometer. In contrast, due to the unavailability of a fast phase modulator at 780 nm, the phase is not randomly modulated but scanned continuously from 0 to $`2\pi `$ using a piezoelectric transducer (PZT). For such a deterministic phase variation, the security of the protocol is not warranted, and thus no genuine secret key could be distributed strictly speaking. However, the experiment provides realistic data, having exactly the awaited structure provided that random phase permutation on Bob’s data are performed.

Due to an imbalance between the paths of the interferometer which modulates the amplitude of the signal beam, the extinction is not strictly zero. In the present experiment that is only aimed at a proof of principle, we subtract the offset field from the data received by Bob. In a real cryptographic transmission, the offset field should be compensated by Alice, either by adding a zeroing field, or by using a better modulator.

All voltages for the electro-optic modulator or the PZT are generated by an acquisition board (National Instruments PCI6111E) connected to a computer. Although all discussions assume the modulation to be continuous, digitised voltages are used in practice. With our experimental parameters, a resolution of 8 bits is enough to hide the amplitude or phase steps under the shot noise. Since the modulation voltage is produced using a 16 bits converter, and the data is digitised over 12 bits, we may fairly assume the modulation and measurement to be continuous.

The homodyne detection was checked to be shot-noise limited for LO power up to $`5\times 10^8`$ photons/pulse. In the present experiment, we used $`1.3\times 10^8`$ photons/pulse for LO power, while each signal pulse contains up to 250 photons. Depending on the run, the overall detection efficiency is either 0.81 or 0.84, due to optical transmission (0.92), mode-matching visibility (0.96 or 0.99) and photodiode quantum efficiency (0.92).

The experiment is thus carried out in such a way that all useful parameters can be measured experimentally. Reconciliation and privacy amplification protocols can thus be performed in realistic – though not fully secret – conditions. The limitations of the present set-up are essentially due to the lack of appropriate fast amplitude and phase modulators at 780 nm. This should be easily solved by operating at telecom wavelengths (1540-1580 nm) where such equipment is readily available. Let us point out also that it is not convenient to transmit separately the signal and LO, so a better solution would be to use a frequency sideband technique similar to Mérolla et al. S1 . Then all light pulses are transmitted together along the same fibre, and a separate radio-frequency is sent from Alice to Bob in order to reconstruct the optical phase information.

### Hypothesis about the detector’s noise : “realistic” vs “paranoid” assumptions

After the quantum exchange, Alice and Bob reveal a subset of their values taken randomly to evaluate the transmission $`G`$ and the total added noise variance. This variance has four contributions: the shot noise $`N_0`$, the channel noise $`\chi _{\mathrm{line}}N_0`$, the electronics noise of Bob’s detector ($`N_{\mathrm{el}}=0.33N_0`$), and the noise due to imperfect homodyne detection efficiency ($`N_{\mathrm{hom}}=0.27N_0`$). The two detection noises $`N_{\mathrm{el}}`$ and $`N_{\mathrm{hom}}`$ originate from Bob’s detection system, so one may reasonably assume that they do not contribute to Eve’s knowledge. This “realistic” assumption has been followed in the article. In that case, the noise from Bob’s detection system also affects Eve’s information so, in equation (4b), we take $`(V_{B|E})_{\mathrm{min}}=N_0/\{G_{\mathrm{line}}(\chi _{\mathrm{line}}+V^1)\}+N_{\mathrm{el}}+N_{\mathrm{hom}}`$, where $`G_{\mathrm{line}}`$ stands for the line transmission.

In contrast, in a “paranoid” approach, one should assume that the noises $`N_{\mathrm{el}}`$ and $`N_{\mathrm{hom}}`$ are also controlled by Eve, that gives her a supplementary advantage. In that case, $`(V_{B|E})_{\mathrm{min}}`$ will be given by $`N_0/\{G(\chi +V^1)\}`$, where $`G`$ now includes both the line and detection efficiencies and $`\chi `$ includes both the line and detection noises. In all cases, the value of $`I_{BA}`$ is given by equation (4a), where $`\chi `$ is the total equivalent noise including both transmission and detection.

Presently we are able to extract practically a key with up to 3.1 dB losses under the “realistic” approach with a reverse reconciliation protocol. Considering now the “paranoid” assumption and reverse reconciliation, the ideal secret key rate is 420 kbits s-1 for a lossless line, and 200 kbits s-1 for $`G_{\mathrm{line}}=0.79`$ (1.0 dB). However, secret bits could be delivered only in the lossless case, at a practical rate of 195 kbits s-1. It is clear that an increase in the reconciliation efficiency would immediately translate into a larger achievable range. Let us point out that we always assume in both the “realistic” or the “paranoid” approach that Eve has an ideal software, quantum memories, perfectly entangled beams, etc. If any of these hypotheses is relaxed, the practical secure range may be extended over the “threshold” presently set by the limited reconciliation efficiency. However, it is not the purpose of the present paper to discuss such “constrained attacks”.

### Implementation of secret key distillation

Secret key distillation was performed by a computer program written in standard C++ that implements the steps described in the paper. Although Alice’s and Bob’s data are both processed on the same computer, it is done in the same way as if the parties were distant and were using a network connection as classical channel. The particular platform used is a regular PC running Linux.

As bursts of data are input to the program, a part of the gaussian key elements are sacrificed and used to estimate the characteristics of the quantum channel. This includes the variances and the correlation coefficient between both sides, which would be exchanged between Alice and Bob over the public authenticated classical channel in a real-life setup.

Depending on the value of the estimated $`I_{AB}`$, the two parties agree on appropriate binary functions (slices 8 ; 10 ) that will transform their gaussian values into bits. These bits are then reconciled, as described in the paper, with sliced error correction 8 ; 10  and an implementation 28  of Cascade 27  as a sub-routine. In our implementation, 5 binary functions are used per gaussian key elements, out of which 2 or 3 (depending on $`I_{AB}`$) are fully disclosed, while the remaining 3 or 2 are reconciled using Cascade.

Next, the data are moved to the privacy amplification routine. Excluding the bits that are fully disclosed and from which no secret key can be extracted, the reconciled bits are processed by use of a transformation randomly taken in a universal class of hash functions 30 ; S7 , which in our case is the class of truncated linear functions in a finite field. First, we consider the reconciled bits as coefficients of a binary polynomial in a representation of the Galois field GF($`2^{110503}`$), hereby called the reconciled polynomial. Then Alice and Bob publicly and randomly choose a random element of the same field and multiply the reconciled polynomial with this chosen element. Finally, they extract from the resulting polynomial the desired number of least significant bits. In our implementation, the representation of the field is GF(2)$`[x]/(p)`$, where $`p=x^{110503}+x^{519}+1`$ is an irreducible polynomial over GF(2), see ref. S8 . The fact that this operation can be implemented efficiently S9  motivated our choice. The size of the field allows us to process up to 110503 bits at once, or equivalently blocks of about 55200 gaussian key elements when Cascade operates on 2 bits per gaussian key element or of 36800 elements with 3 bits per element. To produce a longer key, the gaussian key elements must thus be grouped into blocks.

As explained in the paper, the number of bits that are destroyed by privacy amplification depends on the amount of information that could be inferred by a potential eavesdropper. An easvesdropper Eve has two sources of knowledge. First, she may have attacked the quantum channel and second, she knows the error positions of the reconciled bits from listening to the execution of Cascade. Let $`K`$ be the final key, $`E`$ the ancilla Eve uses for quantum eavesdropping, and $`\mathrm{\Delta }`$ the error positions revealed during reconciliation. We thus need to evaluate $`I(K;E,\mathrm{\Delta })=I(K;E)+I(K;\mathrm{\Delta }|E)`$. The first term on the rhs is upper bounded by $`I_{BE}`$ (in RR) or $`I_{AE}`$ (in DR), while the second term is evaluated numerically for an entangling cloner attack.

This numerical evaluation of $`I(K;\mathrm{\Delta }|E)`$ comes down to integrating $`I(K;\mathrm{\Delta }|E=e)`$ for all possible outcomes $`e`$ of $`E`$, weighted by the probability density function $`p(e)`$ of $`E`$. In the case of the entangling cloner attack, $`E`$ refers both to the knowledge of Eve’s half of the EPR state she injects and to her eavesdropping of the state being sent to Bob, so $`E`$ denotes a bivariate gaussian variable whose covariance matrix can be calculated from the channel characteristics (i.e., attenuation and added noise amplitude). For a given outcome $`e`$ of $`E`$, Eve can infer $`A`$ and $`B`$, Alice’s and Bob’s key elements as a bivariate gaussian variable. Since $`K`$ and $`\mathrm{\Delta }`$ are discrete functions of only $`A`$ and $`B`$, the probability distribution of $`K(A,B)`$ and $`\mathrm{\Delta }(A,B)`$ conditionally on $`E=e`$ can be calculated, hence giving $`I(K;\mathrm{\Delta }|E=e)`$.

Finally, a part of the generated key is used to encrypt 29  the execution of the reconciliation for the next block. The remaining bits, namely the net secret key, can be then used for instance to encrypt the classical communications between Alice and Bob using a one-time pad.


text:
# Evolution of the Colour–Radius Relation and the Morphology–Radius Relation in the SDSS Galaxy Clusters

## 1 Introduction

It is a remarkable feature of galaxy population that there exists a correlation between galaxy type-mix and environment. The correlation between morphology and environment has been studied by many authors (Dressler 1980; Postman & Geller 1984; Whitmore et al. 1993; Whitmore 1995; Dressler et al. 1997; Hashimoto & Oemler 1999; Fasano et al. 2000; Tran et al. 2001; Domínguez et al. 2001, 2002; Helsdon & Ponman 2003; Treu et al. 2003; Goto et al. 2003e). There is also a correlation between colours (spectra) and environment. It is known that the fractions of blue galaxies increases toward outer parts of clusters (e.g., Butcher & Oemler 1984; Rakos & Schombert 1995; Margoniner & de Carvalho 2000; but see Andreon et al. 2003).

Various physical mechanisms have been proposed to explain the correlations. Possible mechanisms include ram-pressure stripping of gas (Gunn & Gott 1972; Farouki & Shapiro 1980; Kent 1981; Fujita & Nagashima 1999; Abadi, Moore & Bower 1999; Quilis, Moore & Bower 2000); galaxy harassment via high speed impulsive encounters (Moore et al. 1996, 1999; Fujita 1998); cluster tidal forces (Byrd & Valtonen 1990; Valluri 1993; Fujita 1998; Gnedin 2003a,b) which distort galaxies as they come close to the centre; interaction/merging of galaxies (Icke 1985; Lavery & Henry 1988; Mamon 1992; Makino & Hut 1997; Bekki 1998; Finoguenov et al. 2003a); evaporation of the cold gas in disc galaxies via heat conduction from the surrounding hot ICM (Cowie & Songaila 1977; Fujita 2003); and a gradual decline in the SFR of a galaxy due to the stripping of halo gas (strangulation or suffocation; Larson, Tinsley & Caldwell 1980; Bekki et al. 2002; Kodama et al. 2001; Finoguenov et al. 2003b). Although these processes are all plausible, however, there exists little evidence demonstrating that any one of these processes is actually responsible for driving galaxy evolution. Since most of these processes act over an extended period of time, observations at a certain redshift cannot easily provide the detailed information that is needed to elucidate subtle and complicated cluster-related processes. It is worth noting that E+A (k+a or post-starburst) galaxies, which often thought to be cluster-related, are found to have their origin in merger/interaction with accompanying galaxies (Goto et al. 2003c,d), and thus E+A galaxies are not likely to be a product of the morphological transition in cluster regions.

To shed light on the underlying physical mechanisms, it is important to trace the redshift evolution of correlations among morphology, colour and environment. Studies on the evolution of morphology–density relation were attempted by Dressler et al. (1997) and Fasano et al. (2000), who reported the significant decrease in S0/E ratio toward increasing redshift up to $`z0.4`$. However, such studies were hampered by the inherent difficulty due to small and inhomogeneous cluster samples (e.g., Andreon 1998). They compared many nearby poor clusters observed with ground-based telescopes with the 10 rich clusters at high redshift observed with the Hubble Space Telescopes. On the colour evolution of cluster galaxies, it has been known that the fractions of blue galaxies in a cluster increases toward higher redshift (the Butcher–Oemler effect; Butcher & Oemler 1978,1984; Rakos & Schombert 1995; Margoniner & de Carvalho 2000; Margoniner et al. 2001; Goto et al. 2003a). To extract more information, it is important to investigate the mutual relationship between the correlations among colour, morphology, and environments.

With the advent of the Sloan Digital Sky Survey (SDSS; Stoughton et al. 2002), we now have a chance to step forward. One of the largest cluster catalogs with a well defined selection function is compiled from the SDSS (Goto et al. 2002a,b). In addition to the availability of the five optical colours in the SDSS ($`u,g,r,i,`$ and $`z`$; Fukugita et al. 1996), this catalog has an advantage in that both high and low redshift clusters are detected using the same cluster finding algorithm (Goto et al. 2002a). In this study, we adopt the cluster-centric-radius as a parameter of environment, and investigate the evolution of colour–radius relation and morphology–radius relation.

The paper is organized as follows: In Section 2, we describe the SDSS data we used. In Section 3, we explain automated morphological classifications and the method to estimate type fraction. In Section 4, we present the results and discuss the physical implications. The cosmological parameters adopted throughout this paper are $`H_0`$=75 km s-1 Mpc-1, and ($`\mathrm{\Omega }_m`$,$`\mathrm{\Omega }_\mathrm{\Lambda }`$,$`\mathrm{\Omega }_k`$)=(0.3,0.7,0.0).

## 2 The SDSS Data

The data used in this study are essentially the same data as used in Goto et al. (2003a). The galaxy catalog is taken from the Sloan Digital Sky Survey Early Data Release (SDSS EDR; Stoughton et al. 2002), which covers $``$394 deg2 of the sky (see Fukugita et al. 1996; Gunn et al. 1998; Lupton et al. 1999,2001,2002; York et al. 2000; Eisenstein et al. 2001; Hogg et al. 2001; Blanton et al. 2003; Richards et al. 2002; Stoughton et al. 2002; Strauss et al. 2002; Smith et al. 2002; Pier et al. 2003 for more detail of the SDSS data). The cluster catalog is presented in Goto et al. (2002a). In this study, we limit our sample clusters to those with richness between 20 and 80, using the cluster richness measured as a number of galaxies between $`M_r^{}=24.0`$ and $`19.44`$ within 0.7 Mpc from a cluster centre after fore/background subtraction (see Goto et al. 2002b). These richness cuts nicely select rich clusters that allows us to study radial profile, and at the same time to reject richest clusters minimizing a possible richness related bias. The resulting number of clusters in our sample is still large ($`N=736`$). We separate these clusters in three redshift bins; $`0.02z0.14`$, $`0.14<z0.20`$, and $`0.20<z0.30`$ (corresponding lookback time of 0.9, 1.9 and 2.7 Gyr). Each redshift bin has 82, 188 and 466 clusters, respectively.

## 3 Analysis

In order to measure blue/late-type fractions, we have to subtract fore/background galaxy counts statistically since our galaxy catalog is based on the imaging data. The adopted procedure is the same to the one used in Goto et al. (2003a) except the use of global galaxy number counts. Instead of local galaxy number counts, we use the entire 394 deg2 of the imaging data to calculate an average galaxy number counts in the magnitude range between $`M_r^{}=24.0`$ and $`20.7`$. The use of global fore/background estimation is essential to probe slight over-densities in cluster perimeter regions. Among the entire region, cluster regions are only 2.4% in angular area. We re-scale this average galaxy number counts to meet an angular area of each radial bin, then subtract fore/background galaxy counts statistically.

Since a physical size of a cluster differs cluster by cluster, it is important to normalize the cluster radius when measuring blue/late type fractions. Since it is difficult to measure a virial radius using the imaging data, we measure a pseudo-virial radius in a following way to normalize the cluster-centric radius (see Goto et al. 2003a for details).

$$R_{vir}=0.7\times (Richness/32)^{1/3}(Mpc)$$
(1)
Equation (1) returns us with a normalized cluster radius of $`R_{vir}=0.7`$ Mpc for clusters with $`Richness=32`$ (c.f., 0.7 Mpc $``$ a virial radius of a system with $`\sigma =`$ 350 km s-1). We note that $`R_{vir}`$ is a pseudo-virial radius, and that caution is needed when comparing with a virial radius from a more realistic measurement such as velocity dispersion or X-ray radial profile. Using a position of the brightest cluster galaxy given in Goto et al. (2002a) as a cluster centre, we sum the number of blue/late and all galaxies within radial bins separated at $`R_{cc}`$/$`R_{vir}=0.2,0.4,0.7`$ and $`1.0`$ ($`R_{cc}`$ is a distance from a cluster centre.), and then subtract fore/background statistically to measure blue/late-type fractions.

To classify galaxies into blue/late-type and others, we use exactly the same method as in Goto et al. (2003a), which we briefly summarize below. We use galaxies between $`M_r^{}=24.0`$ and $`20.7`$. This absolute magnitude range assures the proper measurement of $`ur`$ and $`gr`$ colours for the faintest galaxies at the highest redshift bin. Fractions of blue galaxies are defined in two different ways; $`f_b`$ and $`f_{ur<2.2}`$. The first method, $`f_b`$, uses restframe $`gr`$ colour. If a galaxy is bluer by more than 0.2 in $`gr`$ than the colour of the red-sequence of elliptical galaxies at the redshift of the cluster, we define the galaxy as blue. This definition is a similar definition as was used in original Butcher & Oemler (1984) and recently in Margoniner et al. (2001). The fraction of blue galaxy, $`f_b`$, is defined as the ratio of a number of blue galaxies to that of all galaxies after the global fore/background subtraction. The second method, $`f_{ur<2.2}`$, defines a blue galaxy as a galaxy with $`ur<2.2`$ (Strateva et al. 2001). Again, the fraction of blue galaxy, $`f_{ur<2.2}`$, is defined as the ratio of a number of blue galaxies to that of all galaxies after the global fore/background subtraction.

To separate galaxies into early/late-type morphologically, we use two different purely morphological parameters as are used in Goto et al. (2003a). The first one, $`f_{Cin}`$, uses a concentration parameter, $`Cin`$, defined as the ratio of Petrosian 50% flux radius to Petrosian 90% flux radius. We regard a galaxy with $`Cin>0.4`$ as late-type (Shimasaku et al. 2001). Note that $`f_{Cin}`$ can be affected by the seeing size. Goto et al. (2003a) estimated that an increasing seeing size might increase late-type fractions by $``$5% between $`z=0.3`$ and $`z=0.02`$. The other morphological classification, $`f_{exp}`$, uses galaxy radial profile fitting. When an exponential profile fits better than de Vaucouleurs’ profile to a galaxy, the galaxy is classified into late-type. The fraction of late-type galaxy ($`f_{Cin}`$ or $`f_{exp}`$) is defined as the ratio of a number of late-type galaxies to that of all galaxies after the global fore/background subtraction.

The errors on blue/late-type fractions are computed based on equation (2) of Goto et al. (2003a). The errors are as small as the dot size in Figure 1. Note that these errors do not consider errors in classifying galaxies into sub-categories.

## 4 Results & Discussion

We present the evolution of the colour–radius relation and the morphology–radius relation based on the SDSS data in Figure 1, where blue/late type fractions are plotted against the normalized cluster-centric-radius for the three redshift bins. We use four annulus bins for radial direction (each separated at $`R_{cc}`$/$`R_{vir}=0.2,0.4,0.7`$ and $`1.0`$) and three redshift bins ($`0.02z0.14`$, $`0.14<z0.20`$, and $`0.20<z0.30`$). The $`f_b,f_{ur<2.2},f_{Cin},`$ and $`f_{exp}`$ are plotted in the lower left, upper left, upper right, and lower right panels, respectively; i.e., the two left panels probe colour evolution of cluster galaxies and the two right panels are for morphological evolution. The following observational results can be recognized from Figure 1.

In all the panels, the fractions of blue/late type galaxies decrease with decreasing redshift at almost all cluster-centric-radius. The decrease is consistent with the Butcher–Oemler effect (left panels) and the morphological Butcher–Oemler effect (right panels).

The amount of decrease in fractions might be smaller for colour evolution (left panels; e.g., less than 20% decrease between the highest and the lowest redshift bins) than for morphological evolution (right panels; greater than 30% decrease between the highest and the lowest redshift bins). If the trend is real, it might imply that some red, late-type galaxies changed their morphology into early-type between $`z=0.3`$ and $`z=0.02`$.

Furthermore, at the intermediate redshift ($`0.14<z0.20`$), the fractions of blue galaxies in the two left panels are almost as low as those at the lowest redshift bin ($`0.02z0.14`$). On the other hand, morphological late-type fractions shown in the two right panels are $``$10% larger than those at the lowest redshift bin. The trend might suggest that the colour (spectral) evolution is already completed at the intermediate redshift bin, while the morphological evolution continues to the lowest redshift bin; i.e., timescale of the colour (spectral) evolution may be shorter than the timescale of the morphological evolution.

The radial profile are smooth for the fractions of blue galaxies at any redshift bin. However, for late-type fractions in the two right panels, the radial profiles seem to have a break at $`R_{cc}`$/$`R_{vir}0.3`$ at the highest redshift bin. Then, the profiles become smoother with decreasing redshift. A similar trend is found by Treu et al. (see their Fig. 17) in comparing the morphology–density relation of their cluster Cl0024$`+`$16 with that of the nearby cluster sample in Dressler (1980). The result indicates that the delayed morphological evolution occurs around $`R_{cc}`$/$`R_{vir}=`$ 0.3–0.8, where many spectrally evolved, but morphologically unevolved galaxies possibly exist. The interpretation is consistent with Goto et al. (2003b), who found passive spiral galaxies in the infalling region (perimeter) of clusters. Note that Goto et al. (2003b) used galaxies with spectroscopy whereas this work statistically determines galaxy fractions from the imaging data (see also Couch et al. 1998; Poggianti et al. 1999 for passive spirals).

The difference in timescale of the spectral evolution and the morphological evolution is further clarified in Figure 2, where we plot fractions of four different types of galaxies as a function of cluster-centric-radius. Here, we classified galaxies into the following four different subsamples: blue-late ($`ur<2.2`$, $`Cin0.4`$), blue-early ($`ur<2.2`$, $`Cin0.4`$), red-late ($`ur2.2`$, $`Cin<0.4`$), and red-early ($`ur2.2`$, $`Cin<0.4`$). In the lower right and upper left panels, we plot the redshift evolution of blue-late and red-early types of galaxies as a function of cluster-centric-radius. As expected, red-early type shows gradual increase with decreasing redshift and blue-late type shows gradual decrease with decreasing redshift. However, a contrasting result is found when we compare red-late fractions (upper right) with blue-late fractions (lower right). The fractions of red-late type are not very much dependent on the cluster-centric-radius. In addition, red-late fractions show much evolution between the lowest and the intermediate redshift bins, whereas blue-late fractions show large evolution only between the highest and intermediate redshift bins. Assuming that blue galaxies evolve into red, and that late type evolves into early type, these results suggest that colour evolution (blue to red) mainly occurs between the highest and intermediate redshift, and that morphological evolution (late to early) continues toward the lowest redshift even after colour evolution was completed. In other words, it might indicate that morphological evolution needs longer timescale to be completed ($``$2 Gyr) than the colour evolution does ($``$1 Gyr). These estimates will provide a good benchmark to compare with semi-analytic simulations (e.g., Okamoto & Nagashima 2001; Diaferio et al. 2001; Benson et al. 2001; Springel et al. 2001; Okamoto & Nagashima 2003). Finally, we would like to mention that red-late type galaxies may be the same galaxy population as passive spiral galaxies found in Goto et al. (2003b).

## Acknowledgments

We thank the anonymous referee for many insightful comments which improved the manuscript significantly. T. G. acknowledges financial support from the Japan Society for the Promotion of Science (JSPS) through JSPS Research Fellowships for Young Scientists.