text:
# Equilibrium cluster phases and low-density arrested disordered states: The role of short-range attraction and long-range repulsion

,

## Abstract

We study a model in which particles interact with short-ranged attractive and long-ranged repulsive interactions, in an attempt to model the equilibrium cluster phase recently discovered in sterically stabilized colloidal systems in the presence of depletion interactions. At low packing fraction particles form stable equilibrium clusters which act as building blocks of a cluster fluid. We study the possibility that cluster fluids generate a low-density disordered arrested phase, a gel, via a glass transition driven by the repulsive interaction. In this model the gel formation is formally described with the same physics of the glass formation.

Gels are disordered arrested systems capable of sustaining weak shear stress. Although commonly used in food, cosmetics and detergents, the physical processes involved in the transformation of a solution capable of flowing into a highly viscous substance, and eventually into a gel state, are far from being understood. The gel state can be generated via several routes. Interrupted phase separation, bond percolation, irreversible diffusion limited aggregation, and attractive glass transition have been suggested as possible mechanisms for gel formation in polymeric and colloidal systems.

Sterically stabilized colloidal systems have often been used as an experimental realization of the hard-sphere model, to study both fluid-crystal transition and glass formation. The addition of small polymers in solution modifies the effective colloid-colloid interaction by depletion mechanisms, complementing the hard-sphere repulsion with a short-range attraction, on a length scale of the order of a few percent of the colloid size.

At small volume fractions $`\varphi `$ and low temperature $`T`$, the systems experience dynamic arrest in a transient gel structure poon97 ; lekker-exp-gel . Recent experiments suggest that the gel region in the ($`\varphi T`$) plane is contiguous to a cluster phase, in which particles prlweitz ; daveweitz ; poon-exp-gel  form freely diffusing clusters of finite size. Studies of the phase diagram of short-ranged (uncharged) attractive colloidal systems show that these systems are characterized by fluid and crystal phases and by a very flat metastable liquid-gas spinodal, which extends to considerably low values of $`\varphi `$ lekker ; lekkerkerker02 ; frenkeltermo . In simulations of short-ranged attractive potentials, gel phases have only been observed as a result of an arrested phase separation, but with no contiguity with an equilibrium cluster phase capri .

In this Letter we provide evidence, via molecular dynamics (MD) simulations, that when short-range attractive interactions are complemented by long-range repulsive interactions, equilibrium cluster phases, as well as disordered arrested states contiguous to cluster phases, can be generated. The long-range repulsive interactions may arise from screened electrostatic interactions or from effective interactions accounting for the presence of co-solutes in solution likos ; loewen .

The possibility of generating a cluster phase, i.e., an equilibrium state in which colloidal particles partition into stable clusters via the competition between attraction (which favors cluster growth at low $`T`$) and long-range repulsion (which favors low local particle densities and, therefore, small aggregates) has been discussed in previous theoretical works chandler ; xxxx ; kegel ; muratov . Here we show that in the presence of an appropriate balance between attraction and repulsion the system forms stable clusters, around an optimal cluster size. To quantify this proposition we calculate the ground state energy of clusters of different sizenoi ; wales97  for the case of particles interacting with a potential $`V(r)`$ which is the sum of a short-range attractive potential $`V_{SR}`$ (modeled for simplicity with the generalization to very large $`\alpha `$ of the Lennard-Jones $`2\alpha \alpha `$ potential, as proposed by Vliegenthart et al. lekker )

$$V_{SR}(r)=4ϵ\left[\left(\frac{\sigma }{r}\right)^{2\alpha }\left(\frac{\sigma }{r}\right)^\alpha \right],$$
(1)
and of a screened electrostatic repulsive interaction $`V_Y`$ (modeled by a Yukawa potential)

$$V_Y(r)=A\frac{e^{r/\xi }}{r/\xi }.$$
(2)
In the following, the parameters $`ϵ`$ and $`\sigma `$ entering in $`V_{SR}`$ are chosen as units of energy and distance. The parameter $`\alpha `$ has been chosen equal to $`100`$, corresponding to an attractive range of about three percent of $`\sigma `$ (see inset in Fig.1). For this value of $`\alpha `$, the liquid-gas critical temperature for $`V_{SR}`$ is found to be $`0.235\pm 0.005`$. Fig. 1 shows that the cluster ground state energy, for particles interacting with the potential $`V_{SR}+V_Y`$ for appropriate values of $`A`$ and $`\xi `$ has a minimum at a finite size $`N^{}`$, indicating that clusters of size larger than $`N^{}`$ are energetically disfavored (and, hence, that liquid condensation is inhibited). Since entropic contributions to the free energy will always favor small clusters, aggregation will never proceed beyond $`N^{}`$. The existence of an optimal cluster size (which will of course depend on $`T`$ and $`\varphi `$) suggests that, at low $`T`$, clusters act as building blocks of the supra-molecular cluster fluid.

The cluster-cluster interaction is essentially controlled by the Yukawa repulsion, due to the short-range nature of the attractive part of the potential. It can be calculated analytically, if one approximates the cluster to a sphere of radius $`R`$ composed of a homogeneous density $`\rho `$ of points interacting via the potential $`V_Y`$ (Eq. (2)). The resulting cluster-cluster interaction is, again, of the Yukawa type with the same $`\xi `$ but a renormalized amplitude $`A(R)`$

$$\frac{A(R)}{A}=\left\{2\pi \xi ^3\rho e^{R/\xi }\left[1+\frac{R}{\xi }+\left(\frac{R}{\xi }1\right)e^{2R/\xi }\right]\right\}^2.$$
(3)
Therefore, aggregation of particles into clusters of size $`N`$ (and radius $`RN^{1/3}`$) has two effects: i) renormalization of the amplitude of the interaction, and ii) decrease of the cluster number density $`n`$ according to

$$\frac{n(R)}{n}=\frac{1}{N}.$$
(4)
Since $`\xi `$ is not modified, the aggregation process can be seen as the generator of a flow, defined via Eqs. (3) and (4), in the phase diagram of the pure repulsive Yukawa potential.

Fig. 2 shows an extended phase diagram for the pure Yukawa potential. Together with the fluid-crystal coexistence lines, previously calculated by Kremer et al. grest , and Meijer and Frenkel frenkel , we show an iso-diffusivity curvesimyuk  (i.e., the locus of point with the same diffusion constant) and the ideal glass line, evaluated within the Mode Coupling Theory mct  (MCT) using exact numerical structure factors as input. The latter curve provides an upper limit to the glass transition line. It is important to note that the long-range nature of the repulsion enables the formation –at low $`k_BT/A`$ values– of solid phases with extremely low density (named Wigner crystal or Wigner glass bosse ; lai97 ; laponite ).

An example of the flow in the Yukawa phase diagram induced by clustering is shown in Fig. 2. The starting point is assumed to be a system composed of monomers ($`N=1`$) interacting via the potential $`V_{SR}+V_Y`$. If $`T`$ is sufficiently low, a clustering phenomenon will take place, driven by the attractive part of the potential. Such aggregation process acts as a generator of a flow in the phase diagram, moving the system, due to the effective increase of the repulsive amplitude (Eq. 3) and to the associated decrease of the cluster number density (Eq. 4), closer and closer to the melting and glass lines. If the optimal size is reached before crossing the melting line, the equilibrium state for the system will correspond to a fluid phase of clusters. In the other case, an arrest phenomenon will take place, driven by the crystallization or –more likely because of the inherent self-generated cluster polydispersity– by a glass transition. With this mechanism arrested disordered states with small $`\varphi `$ can be generated. This new route to the gel state demonstrates, in the case of weakly charged short-ranged attractive colloids, the identity of the gel and the glass states of matter.

We next provide numerical evidence that the above proposed mechanism does lead to disordered arrested states contiguous to a cluster phase. We perform MD simulations of particles interacting with the potential $`V_{SR}+V_Y`$simudetails , both along a constant $`T`$ and along a constant $`\varphi `$ path. Here, we focus on the specific choice of $`A=0.05`$ and $`\xi =2`$, referring to future publications for a more extensive analysis. This choice of parameters is characteristic of colloidal particles in weakly polar solvents daveweitz ; bladereen ; grier .

The cluster size distribution, shown in Fig. 3, provides a characterization of the clustering process. While at high $`T`$ the cluster size is monotonically decreasing, around $`k_BT/ϵ=0.2`$ the distributions develop a peak at a finite value, signaling the presence of stable equilibrium aggregates composed of $`2030`$ particles. We note that such resulting cluster phase can be considered as an echo of the interrupted spinodal decomposition, since clustering takes places just below the liquid-gas critical temperature of the bare $`V_{SR}`$ potential.

For each of the studied $`T`$, the unclustered system ($`N=1`$) can be represented as a point in Fig. 2. A flow curve, starting from each of these points, can be calculated according to Eqs. 3 and 4. In this way, a theoretical prediction for the average cluster size at which crystallization and (MCT) vitrification are expected can be derived. Such curves are shown in the inset together with the $`T`$-dependence of the average cluster size calculated from the MD data. For the present choice of $`A`$, $`\xi `$ and $`\varphi `$, structural arrest (either in the form of crystallization or in the form of vitrification) is predicted when the average cluster size is $`35`$ particles ($`k_BT/ϵ0.15`$).

Before turning to the evaluation of dynamical properties, we discuss the evolution of the structure with $`T`$, analyzing the static structure factor $`S(q)`$. As shown in Figure 4-a, $`S(q)`$ progressively develops a peak at wave-vectors $`q`$ associated with distances in real space of the order of several particle diameters, reflecting the nearest neighbor cluster-cluster distance. The amplitude of this novel peak changes significantly in a small temperature range around the $`T`$-region where the flow diagram predicts the crossing of the crystallization line. Therefore, the system changes from a solution of interacting particles toward a system of larger and larger interacting clusters. In the cluster phase, the location of the peak changes very weakly both with $`T`$ and with $`\varphi `$. It is interesting to note that both the location and the amplitude of the cluster-cluster peak are in agreement with the experimental data for colloidal gels reported in Ref. prlweitz . To confirm that the aggregation process indeed generates a disordered arrested state, we show in Fig. 4-(Bottom) the coherent scattering function $`F(q,t)`$. It displays a (one-step) decay, with a characteristic time which increases on cooling. At the lowest studied $`T`$, the correlation function does not decay to zero, signaling that the system is trapped in a non-ergodic state. The transition from fluid to glass is very sharp, again in full agreement with results for colloidal gels prlweitz .

In summary, we have proposed and studied a model which shows equilibrium cluster phases and arrested states at low densities, opening the way for accurate comparisons between experiments and theoretical predictions. Contrary to what has been proposed before for attractive colloids, we suggest that, in these systems, the gel phase is not stabilized by the short-range attraction (and hence its formation is not connected to the attractive glass transition observed at high packing fractions fuchs ). The simulation data (Figs. 3 and 4) show that the competition between isotropic short-range attraction and long-range repulsion — of the type found in colloidal systems and in protein solutions piazza00 ; pellicane03  — is capable of generating a very rich range of phenomena douglas . Within the present model, the formation of low density arrested states can be modeled as a glass transition, where clusters (as opposed to particles) are trapped in cages generated by the long-range repulsions. It will be interesting to study if the arrested disordered state is metastable with respect to a a fcc or bcc crystal of ”spherical” clusters (a Wigner crystal) or with respect to lamellar phases, which has been observed in related models chandler ; xxxx ; muratov .

We thank D. Chandler, D. Frenkel, W. Kegel, D.R. Reichman, G. Ruocco, and I. Saika-Voivod for comments. F. S. thanks W. Poon for suggesting to investigate the role of long-range repulsive interactions in short-range attractive colloidal systems. Authors acknowledge support from MIUR Cofin 2002 and Firb.


format: text

text:
# On the Boundary Entropy of One-dimensional Quantum Systems at Low Temperature

## I Introduction

The logarithm of the partition function $`\mathrm{ln}Z=\mathrm{ln}\mathrm{Tr}\mathrm{exp}(\beta H)`$ for a one-dimensional quantum critical system with a boundary takes the universal formAffleck and Ludwig (1991) $`(c\pi /6)(L/\beta )+\mathrm{ln}g`$ where $`H`$ is the hamiltonian, $`\beta =1/T`$ the inverse temperature, $`L\beta `$ is the length, $`c`$ is the numerical coefficient of the bulk conformal anomaly, and $`g`$ is the¡ “universal noninteger ground-state degeneracy” at the boundary (using natural units in which $`\mathrm{}=k=v=1`$, $`v`$ being the velocity of “light”). This formula applies in the limit of large $`L`$. The number $`g`$ is an invariant of the universality class of the critical boundary condition. It was conjectured that $`g`$ decreases from critical point to critical point under renormalizationAffleck and Ludwig (1991, 1993).

For a 1-d quantum system that is critical in the bulk but is *not* critical at the boundary, the logarithm of the partition function at low temperature can be written in the form $`\mathrm{ln}Z=(c\pi /6)(L/\beta )+\mathrm{ln}z_L`$ and the boundary partion function $`z`$ can be defined as $`lim_L\mathrm{}z_L`$. That is, the partition function takes the universal form

$$\mathrm{ln}Z=(c\pi /6)(L/\beta )+\mathrm{ln}z$$
(1)
up to corrections which vanish in the limit $`L\mathrm{}`$.

Quantum critical points occur at zero temperature. At temperature $`T>0`$, the correlation functions are of course not scale-invariant, but decay exponentially at distance scale $`\beta =1/T`$. Nevertheless, it is meaningful to talk of a bulk critical system at temperature $`T>0`$. The low energy degrees of freedom and their coupling constants can be identified in the $`T=0`$ critical system, which is described by a 1+1 dimensional quantum field theory. That quantum field theory, at temperature $`T>0`$, describes the quantum critical system at temperature $`T`$. It is possible to identify temperature-independent coupling constants within the critical system. These coupling constants can be held to their critical values when the temperature is above zero. Quantum critical phenomena are special in this respect. In generic critical phenomena, which have non-zero critical temperatures, the temperature is just one of the coupling constants.

Given that the bulk system is critical and that $`L=\mathrm{}`$, the only dimensionful parameter is the temperature $`T`$. The logarithm of the boundary partition function is thus a function $`\mathrm{ln}z(\mu \beta )`$ that depends only on the temperature, in units of $`\mu `$, where $`\mu `$ is a small temperature that sets the renormalization scale (or, equivalently, a small energy or inverse time or inverse distance). The total entropy then takes the universal form

$$S=(1\beta /\beta )\mathrm{ln}Z=(c\pi /3)(L/\beta )+s(\mu \beta )$$
(2)
where $`s(\mu \beta )=(1\beta /\beta )\mathrm{ln}z`$ is the boundary entropy. At a critical point, $`s`$ is equal to the constant $`\mathrm{ln}g`$.

We prove here a gradient formula

$$g_{ab}(\lambda )\beta ^b(\lambda )=s/\lambda ^a$$
(3)
where the $`\lambda ^a`$ form a complete set of boundary coupling constants, $`g_{ab}(\lambda )`$ is a certain metric on the space of all boundary conditions, and $`\beta ^a(\lambda )`$ is the boundary $`\beta `$-function. It follows directly from the gradient formula that $`\mu s/\mu =\beta ^a_as=g_{ab}\beta ^a\beta ^b`$ so $`s`$ decreases under the renormalization group except where $`\beta ^a=0`$, at the critical points. The gradient formula eliminates the possibility of esoteric asymptotic behavior under renormalization. Recurring trajectories such as limit cycles are excluded. The $`g`$-conjecture for the rg flows between critical points is a corollary of the gradient formula.

The gradient formula implies equally that the boundary entropy decreases with temperature, $`\beta s/\beta =\mu s/\mu <0`$. The total entropy $`S`$ obviously decreases with temperature, because $`S/\beta =\beta (HH)^2`$. However, the decrease of the bulk contribution $`(\pi /3)(L/\beta )c`$ masks the change in $`s`$, so it is not obvious that the boundary entropy by itself must decrease with temperature. The gradient formula implies that it does. It follows that the thermodynamic boundary energy also decreases with temperature, $`\mu \beta u/\beta =s/\beta <0`$.

Complete control over the possible behavior at asymptotically low temperature is still lacking, because we do not prove that $`s`$ is bounded below. If $`s`$ is bounded below, then the system must go to a critical point at zero temperature. Of course, the total entropy $`S`$ of any system is bounded below, as long as the system is of finite size. So, for any finite size $`L`$, $`s_L=S(c\pi /3)(L/\beta )`$ is bounded below, as $`T0`$. However, the lower bound can descend without limit as $`L\mathrm{}`$, so $`s=lim_L\mathrm{}s_L`$ is not necessarily bounded below, as $`T0`$. It still remains to be proved that $`s(\mu \beta )`$ is bounded below, as $`T0`$. If the boundary entropy is bounded below, then the boundary energy is also bounded below.

The gradient formula that we prove is mathematically equivalent to a gradient formula conjectured in string theoryWitten (1992, 1993); Shatashvili (1993a, b); Kutasov et al. (2000). Evidence was given for the string theory conjectureWitten (1993); Shatashvili (1993b); Konechny (2004), but the formula was never proved. It has been claimed that a proof was given in Ref. Kutasov et al. (2000), but it was assumed there that the boundary $`\beta `$-function $`\beta ^a(\lambda )`$ is linear in the coupling constants $`\lambda ^a`$. This is an invalid assumption. The $`\beta `$-function cannot be linearized when there are marginally relevant couplings or, more generally, whenever resonance conditions occur (as discussed, for example, in Ref. Shatashvili (1993b)). Moreover, the conjectured string gradient formula is expressed in un-physical quantities, in terms of un-normalized correlation functions. Our contribution is to express the gradient formula in terms of normalized correlation functions and the boundary entropy, which are physical quantities, and to prove the formula using physical properties of the 1-d quantum system. Some of the ideas used in the proof can be found in the string theory workWitten (1992); Shatashvili (1993b). The re-writing of the conjectured string gradient formula is based on an idea that is implicit in Ref. Kutasov et al. (2000) and was mentioned explicitly to usMoore . Here, to avoid distracting from the physical meaning, we first prove the gradient formula in physical terms, and only afterwards explain the connection to the string conjecture.

The proof of the gradient formula applies to *all* local 1d quantum systems. It uses only the basic principles of quantum mechanics and locality. The gradient formula must therefore hold in every local 1d quantum mechanical model. The point of proving a result such as the gradient formula is to give reliable theoretical information about what is physically possible. For instance, when building devices out of low temperature 1-d quantum systems joined at boundaries, it will be useful to know in advance, with certainty, what kinds of boundary behaviors are possible. It will be useful to know that the boundary must always behave as a thermodynamic system, except that it does not obey the third law. Proof also reveals what must be done to evade the theoretical limits. The gradient formula itself is not likely to be avoidable, since the proof depends only on the basic principles of quantum mechanics and renormalization, assuming only the existence of a local stress-energy tensor, which is assured by microscopic locality. Rather, attention is directed towards exotic systems, where the metric $`g_{ab}(\lambda )`$ degenerates, or where $`s`$ is infiniteFriedan (1999, 2003); Tseng (2002), or even where $`s`$ might not be bounded below, if this cannot be proved impossible. A lower bound on $`s`$ would have to depend on the details of the bulk system. The bound could not be uniform, not a function of $`c`$ alone. This can be seen in the critical $`c=1`$ gaussian model, where the values of $`g`$ depend on the marginal coupling constant of the bulk model, and can become arbitrarily close to zeroElitzur et al. (1999).

## II The stress-energy tensor in the presence of a boundary

The equilibrium observables of the system live on the cylindrical euclidean spacetime, periodic in euclidean time with period $`\beta `$ (see Fig. 1). The spacetime coordinates are $`x^\mu =(x,\tau )`$, $`0x<L`$, $`\tau \tau +\beta `$. The boundary is at $`x=0`$.

The stress-energy tensor $`T_{\mu \nu }(x,\tau )`$ expresses the response of the system to an infinitesimal local variation of the metric, $`g_{\mu \nu }\delta _{\mu \nu }+\delta g_{\mu \nu }(x,\tau )`$,

$$\delta \mathrm{ln}Z=(1/2)\mu ^2d\tau dx\delta g_{\mu \nu }T^{\mu \nu }(x,\tau ).$$
(4)
We specialize to 1+1 dimensions the general analysis of the stress-energy tensor in space-times with boundaryMcAvity and Osborn (1993). The stress-energy tensor can be written as a bulk part plus a boundary part

$$\delta \mathrm{ln}Z=(1/2)\mu ^2d\tau dx\delta g_{\mu \nu }T_{bulk}^{\mu \nu }(x,\tau )+(1/2)_0^\beta \mu d\tau \delta g_{\mu \nu }(0,\tau )t^{\mu \nu }(\tau )$$

$$T_{\mu \nu }=T_{\mu \nu }^{bulk}(x,\tau )+\delta (\mu x)t_{\mu \nu }(\tau ).$$
(5)
There could also be a boundary operator proportional to $`\delta ^{}(x)\mathrm{𝟏}`$, but the identity operator makes no contribution to connected correlation functions, so we can ignore it.

The conservation equations follow from invariance of the physics under localized coordinate reparametrizations $`\delta x^\mu =v^\mu (x,\tau )`$ where the vector field $`v^\mu `$ is tangent to the boundary, i.e. $`v^x(0,\tau )=0`$. The coordinate reparametrization is equivalent to a change in the metric tensor $`\delta g_{\mu \nu }=_\mu v_\nu +_\nu v_\mu .`$ Plugging this into the formula for $`\delta \mathrm{ln}Z`$ and setting the variation to zero we obtain, after integration by parts, the bulk conservation equation $`^\mu T_{\mu \nu }^{bulk}=0`$ and also

$$\mu d\tau \left(\mu T_{x\tau }^{bulk}v^\tau t_{\mu \nu }^\mu v^\nu \right)=0.$$

at the boundary, which is equivalent to the boundary conservation equations $`t_{xx}=t_{x\tau }=t_{\tau x}=0`$ and

$$\mu T_{x\tau }^{bulk}(0,\tau )+_\tau \theta (\tau )=0$$
(6)
where $`\theta (\tau )t_{\tau \tau }(\tau )`$. The boundary operator $`\theta `$ was described in Ref. Ghoshal and Zamolodchikov (1994).

The trace of the stress-tensor is

$$T_\mu ^\mu =\mathrm{\Theta }(x,\tau )=\mathrm{\Theta }_{bulk}(x,\tau )+\delta (\mu x)\theta (\tau ).$$
(7)
The system is critical in the bulk, so $`\mathrm{\Theta }_{bulk}(x,\tau )=0`$ up to contact terms. The full trace is $`\mathrm{\Theta }=\delta (\mu x)\theta (\tau )`$, entirely a boundary operator.

The space of boundary conditions is parameterized by the coupling constants $`\lambda ^a`$ which couple to the renormalized local boundary fields $`\varphi _a`$

$$_a\mathrm{ln}Z=\mu d\tau \varphi _a(\tau ).$$

The boundary trace $`\theta `$ can be decomposed into a linear combination of the boundary fields and the identity operator

$$\theta =\beta ^a(\lambda )\varphi _a+h(\lambda )\mathrm{𝟏}$$
(8)
where the coefficients $`\beta ^a(\lambda )`$ are the boundary $`\beta `$-functions. We will not have to worry about the term $`h(\lambda )\mathrm{𝟏}`$, because $`\theta `$ will only appear within connected correlation functions.

The foregoing are operator statements. In correlation functions, the stress-energy tensor will also have contact terms. The generator of dilatations $`\delta g_{\mu \nu }=2(\delta \mu /\mu )\delta _{\mu \nu }`$ is $`\mu ^2d\tau dx\mathrm{\Theta }(x,\tau )`$ so the renormalization group equation for $`\mathrm{ln}Z`$ is

$`(\mu /\mu )\mathrm{ln}Z`$ $`=`$ $`{\displaystyle \mu ^2d\tau dx\mathrm{\Theta }(x,\tau )}`$ (9)
$`=`$ $`{\displaystyle \mu ^2d\tau dx\delta (x)\theta (\tau )+\mathrm{\Theta }_{bulk}(x,\tau )_c}`$
$`=`$ $`\beta ^a_a\mathrm{ln}Z+\mu \beta h(\lambda ).`$
For the one-point functions,

$`(\mu /\mu )\varphi _b(\tau _1)`$ $`=`$ $`{\displaystyle \mu ^2d\tau dx\mathrm{\Theta }(x,\tau )\varphi _b(\tau _1)_c}`$ (10)
$`=`$ $`{\displaystyle \mu ^2d\tau dx\left[\delta (x)\theta (\tau )+\mathrm{\Theta }_{bulk}(x,\tau )\right]\varphi _b(\tau _1)_c}`$
$`=`$ $`\beta ^a_a\varphi _b(\tau _1)+(\gamma _b^a\delta _b^a)\varphi _a(\tau _1)`$
where the coefficients $`\gamma _b^a\delta _b^a`$ come from contact terms of $`\mathrm{\Theta }_{bulk}`$ and $`\theta `$ with $`\varphi _a`$. Because of the contact terms, $`\mathrm{\Theta }_{bulk}`$ cannot be omitted. The identity $`\gamma _b^a=_b\beta ^a`$ follows from $`[\mu /\mu ,_a]=0`$, which in turn follows from the definition of the $`\lambda ^a`$ as the coupling constants renormalized at scale $`\mu `$.

We will need one last property of the stress-energy tensor, that $`T_{\mu \nu }^{bulk}(x,\tau )`$ decays as $`\mathrm{exp}(4\pi x/\beta )`$ in connected correlation functions far from the boundary. When $`x`$ is far from the boundary, $`T_{\mu \nu }^{bulk}(x,\tau )`$ behaves as in the bulk theory without boundary. The exponential decay condition expresses the conformal invariance of the bulk critical system. It is derived using the interpretation of Fig. 1 in which space and imaginary time are exchanged. Space becomes a circle of length $`\beta `$. The correlation functions become the expectation values $`B|\mathrm{}|0`$ where $`|0`$ is the ground state of the bulk critical system on the circle, and $`B|`$ is the state representing the boundary condition. Using the complex coordinate $`w=2\pi (x+i\tau )/\beta `$, the bulk stress-energy tensor takes the form $`T_{xx}^{bulk}=T_{\tau \tau }^{bulk}=T(w)+\overline{T}(\overline{w})`$, $`T_{x\tau }^{bulk}=i(T(w)\overline{T}(\overline{w}))`$, where $`T(w)=_nL_n\mathrm{exp}(nw)`$ and $`\overline{T}(\overline{w})=_n\overline{L}_n\mathrm{exp}(n\overline{w})`$, $`n`$ ranging from $`\mathrm{}`$ to $`+\mathrm{}`$, the $`L_n`$ and $`\overline{L}_n`$ being the Virasoro operators. Bulk conformal invariance means $`L_1|0=\overline{L}_1|0=0`$. Therefore, $`T(w)\mathrm{exp}(2w)`$ and $`T(\overline{w})\mathrm{exp}(2\overline{w})`$ in connected correlation functions, far from the boundary. So $`T_{\mu \nu }^{bulk}(x,\tau )\mathrm{exp}(4\pi x/\beta )`$.

## III The proof

We prove the gradient formula, Eq. 3, with the metric on the space of boundary conditions given by

$$g_{ab}(\lambda )=\mu d\tau _1\mu d\tau \varphi _a(\tau _1)\varphi _b(\tau )_c\left(1\mathrm{cos}\left[2\pi (\tau \tau _1)/\beta \right]\right).$$
(11)
This is essentially the metric proposed in Ref. Witten (1992), except that Ref. Witten (1992) used the un-normalized, full two-point function, while we use the normalized, connected two-point function. Because we are using the connected two-point function, we can write

$$g_{ab}\beta ^b=\mu d\tau _1\mu d\tau \varphi _a(\tau _1)\theta (\tau )_c\left(1\mathrm{cos}\left[2\pi (\tau \tau _1)/\beta \right]\right).$$
(12)
The identity component of $`\theta `$ makes no contribution to the connected two-point function. Let us deal with the term containing the cosine:

$`A_a(\tau _1)`$ $``$ $`{\displaystyle \mu d\tau \varphi _a(\tau _1)\theta (\tau )_c\left(\mathrm{cos}\left[2\pi (\tau \tau _1)/\beta \right]\right)}`$ (13)
$`=`$ $`{\displaystyle \mu d\tau \varphi _a(\tau _1)_\tau \theta (\tau )_c(\beta /2\pi )\mathrm{sin}\left[2\pi (\tau \tau _1)/\beta \right]}`$
We have integrated by parts on the boundary to obtain the second equation. The correlation functions are distributions, so integration by parts is justified. By the boundary conservation law Eq. 6,

$$A_a(\tau _1)=\mu d\tau \varphi _a(\tau _1)\mu T_{x\tau }^{bulk}(\tau )_c(2)v^\tau (0,\tau )$$
(14)
where we define $`v^\tau (0,\tau )(\beta /4\pi )\mathrm{sin}\left[2\pi (\tau \tau _1)/\beta \right]`$ as a tangent vector field on the boundary. Next, we extend the boundary vector field $`v^\tau (0,\tau )`$ to a conformal Killing vector field $`v^\mu (x,\tau )`$ in the bulk. That is, $`v^x(0,\tau )=0`$ and $`_\mu v_\nu +_\nu v_\mu =g_{\mu \nu }_\sigma v^\sigma `$. Such a vector field is most easily found as an analytic vector field $`v^w=(2\pi /\beta )(v^x+iv^\tau )`$ in the complex coordinate $`w=2\pi (x+i\tau )/\beta `$,

$`v^w`$ $`=`$ $`(2\pi /\beta )(v^x+iv^\tau )=\left[\mathrm{exp}(ww_1)\mathrm{exp}(w+w_1)\right]/4`$
$`\overline{v}^{\overline{w}}`$ $`=`$ $`(2\pi /\beta )(v^xiv^\tau )=\left[\mathrm{exp}(\overline{w}\overline{w}_1)\mathrm{exp}(\overline{w}+\overline{w}_1)\right]/4.`$
Then

$$_\sigma v^\sigma =_wv^w+_{\overline{w}}\overline{v}^{\overline{w}}=\mathrm{cos}\left[2\pi (\tau \tau _1)/\beta \right]\mathrm{cosh}(2\pi x/\beta ).$$

Now we integrate by parts in the bulk, using the bulk conservation equation, to obtain

$$A_a(\tau _1)=\mu ^2d\tau dx\varphi _a(\tau _1)T_{\mu \nu }^{bulk}(x,\tau )_c(^\mu v^\nu +^\nu v^\mu ).$$
(15)
There is no boundary term at large $`x`$ because of the decay condition $`T_{\mu \nu }^{bulk}(x,\tau )\mathrm{exp}(4\pi x/\beta )`$. Then we use the fact that $`v^\mu `$ is a conformal Killing vector to write

$$A_a(\tau _1)=\mu ^2d\tau dx\varphi _a(\tau _1)\mathrm{\Theta }_{bulk}(x,\tau )_c_\sigma v^\sigma .$$
(16)
Finally, we can approximate $`_\sigma v^\sigma 1`$, because $`\mathrm{\Theta }_{bulk}=0`$ except for contact terms. The error term is

$$\mu ^2d\tau dx\varphi _a(\tau _1)\mathrm{\Theta }_{bulk}(x,\tau )_c(_\sigma v^\sigma 1).$$

The boundary operator $`\varphi _a(\tau _1)`$ is renormalizable, and $`\mathrm{\Theta }_{bulk}`$ has dimension $`2`$, so the most singular contact terms in the two-point function are of the form $`\delta (x)\delta ^{}(\tau \tau _1)`$ and $`\delta ^{}(x)\delta (\tau \tau _1)`$. But $`_\sigma v^\sigma (x,\tau )1`$ vanishes to second order at $`x=0`$, $`\tau =\tau _1`$, so there is no error. Thus

$$A_a(\tau _1)=\mu ^2d\tau dx\varphi _a(\tau _1)\mathrm{\Theta }_{bulk}(x,\tau )_c.$$
(17)
Using Eq. 17 in Eq. 12, we arrive at

$`g_{ab}\beta ^b`$ $`=`$ $`{\displaystyle \mu d\tau _1\mu ^2d\tau dx\varphi _a(\tau _1)\left[\delta (\mu x)\theta (\tau )+\mathrm{\Theta }_{bulk}(x,\tau )\right]_c}`$ (18)
$`=`$ $`{\displaystyle \mu d\tau _1\mu ^2d\tau dx\varphi _a(\tau _1)\mathrm{\Theta }(x,\tau )_c}`$
$`=`$ $`{\displaystyle \mu d\tau _1(\mu /\mu )\varphi _a(\tau _1)}`$
$`=`$ $`(\mu /\mu 1)_a\mathrm{ln}Z`$
$`=`$ $`_a(\mu /\mu 1)\mathrm{ln}Z`$
$`=`$ $`_as`$
which is the gradient formula.

## IV Comments on the gradient formula

Each element of the gradient formula is covariant under renormalization. The boundary entropy $`s`$ is covariant under renormalization, $`\mu s/\mu =\beta ^a_as`$, even though the partition function is not (see Eq. 9). Using Eq. 9,

$`(\mu /\mu \beta ^a_a)s`$ $`=`$ $`(\mu /\mu \beta ^a_a)(1\mu /\mu )\mathrm{ln}Z`$ (19)
$`=`$ $`(1\mu /\mu )(\mu /\mu \beta ^a_a)\mathrm{ln}Z`$
$`=`$ $`(1\mu /\mu )(\mu \beta h)=0.`$
That is, the entropy is not sensitive to a shift of the ground state energy. The covariance of $`\beta ^a`$ is just its $`\mu `$-independence. The metric $`g_{ab}`$ is covariant under renormalization because it is defined in terms of normalized, connected correlation functions, in Eq. 11.

To show that the metric $`g_{ab}`$ is positively definite, we need only remark that $`g_{ab}\delta \lambda ^a\delta \lambda ^b`$ is given in Eq. 11 as a positive two point function of $`\varphi =\delta \lambda ^a\varphi _a`$, integrated against a positive function.

The cosine term in the metric plays a twofold role. On the one hand, it provides the $`\mathrm{\Theta }_{bulk}`$ term in the correlation function of $`\mathrm{\Theta }`$ with the boundary operator. On the other hand the cosine term renders the metric independent of contact terms in the two-point functions of the boundary operators. Such terms could spoil the positivity of the metric. The metric, as defined by Eq. 11, is independent of contact terms. During the proof of the gradient formula, we split it into two parts, each of which does depend on the contact terms. At that point, the two point functions have to be treated as distributions. In the end, when the two terms are joined together, the result is independent of the contact terms. The technical roles of the cosine term are evident, but we do not see a deeper meaning. The cosine first appeared in the string theory metric proposed in Ref. Witten (1992). But the proposal was not natural in string theory, as it involved integrating dimension zero fields. So we still do not see a natural interpretation of the cosine term.

## V Relation with string theory

The conjectured string theory gradient formula involves an additional boundary coupling constant $`\lambda ^0`$ which couples to the identity operator $`\varphi _0=\mathrm{𝟏}`$. The string partition function is

$$Z_s=\mathrm{exp}(\mu \beta \lambda ^0)z(\mu \beta )$$

where $`z(\mu \beta )`$ is the boundary partition function, from Eq. 1. The string $`\beta `$-function, $`\beta _s^a`$, is the ordinary $`\beta ^a`$ for the ordinary coupling constants, plus, from Eq. 9, $`\beta _s^0=\lambda ^0+h(\lambda )`$. The conjectured string theory gradient formula is

$$G_{ab}^s\beta _s^b=_ag_s$$

where

$$g_s=(1\mu /\mu )Z_s$$

and the string metric is

$$G_{ab}^s(\lambda )=\mu d\tau _1\mu d\tau Z_s\varphi _a(\tau _1)\varphi _b(\tau )\left(1\mathrm{cos}\left[2\pi (\tau \tau _1)/\beta \right]\right).$$
(20)
These string formulas are un-physical, when applied to 1-d quantum systems. No physical probe couples to the identity operator $`\varphi _0=\mathrm{𝟏}`$, so $`\lambda ^0`$ is not a physical coupling constant. Un-normalized correlation functions are not measurable. Changes in $`g_s`$ are not locally measurable, because $`g_s`$ is constructed from $`z`$, not $`\mathrm{ln}z`$. On the other hand, all of the elements of the physical gradient formula, Eq. 3, can by measured by local operations at the boundary of the 1-d system. The string gradient formula *is* formally sensible from the string theory perspective. The $`\lambda ^a`$ are the wave-modes of spacetime fields, $`\lambda ^0`$ is the zero-mode of the tachyon field. The equation $`\beta _s^a=0`$ has the form of a space-time equation of motion. The function $`g_s(\lambda )`$ has the form of a space-time action.

The un-physical parameter $`\lambda ^0`$ can be eliminated by extremizing $`g_s`$Kutasov et al. (2000); Moore . We carry out this idea. We calculate that $`g_s/\lambda ^0=0`$ at $`\lambda ^0=a_0`$, $`\mu \beta a_0=\mu \mathrm{ln}z/\mu `$. We calculate that, at $`\lambda ^0=a_0`$, $`g_s=Z_s=z\mathrm{exp}(\mu \beta a_0)`$, which is the physical quantity $`\mathrm{exp}(s)`$. It now becomes straightforward to show the equivalence between the string gradient formula and the physical gradient formula. The string gradient formula is trivial in the direction of $`\lambda ^0`$, and is precisely the physical formula on the subspace $`\lambda ^0=a_0`$. To be explicit, the components of the string metric are $`G_{00}^s=(\mu \beta )^2Z_s`$, $`G_{0b}^s=\mu \beta Z_s_a\mathrm{ln}Z_s`$, $`G_{ab}^s=Z_s(g_{ab}+_a\mathrm{ln}Z_s_b\mathrm{ln}Z_s)`$, where the indices $`a,b`$ now range only over the physical coupling constants. The string gradient formula splits into two equations

$`G_{00}^s\beta _s^0+G_{0b}^s\beta ^b=_0g_s`$
$`G_{a0}^s\beta _s^0+G_{ab}^s\beta ^b=_ag_s.`$
The first is satisfied identically, it is just the rg equation for $`Z_s`$,

$$\mu Z_s/\mu =(\beta ^a_a+\mu \beta \beta _s^0)Z_s.$$

The second equation, after substituting and then using the rg equation for $`Z_s`$, becomes

$$Z_sg_{ab}\beta ^b=Z_s(1+\mu /\mu )_a\mathrm{ln}Z_s$$

which is exactly the physical gradient formula, since $`_a\mathrm{ln}Z_s=_a\mathrm{ln}Z`$. So, by proving the physical gradient formula, we have also proved the string gradient formula.

We would like to thank Gregory Moore and Alexander Zamolodchikov for stimulating discussions. We especially thank Gregory Moore for pointing out a deficiency in an early version of the proof. This work was supported by the Rutgers New High Energy Theory Center, which A. K. thanks for warm hospitality. The work of A. K. was supported in part by BSF-American-Israel Bi-National Science Foundation, the Israel Academy of Sciences and Humanities-Centers of Excellence Program, the German-Israel Bi-National Science Foundation.


text:
# Preparing Red-Green-Blue (RGB) Images from CCD Data

## 1 Introduction

The purest and most thoughtful minds are those which love colour the most 111John Ruskin, The Stones of Venice

One of the reasons that our forbears were drawn to astronomy was the sheer drama and beauty of the night sky. When photography became a standard astronomical technique, it was natural to combine images taken in different bands to produce spectacular colored images of the night sky (e.g. Malin, (1992) and references therein).

When CCD detectors were introduced into astronomy 25 years ago, and especially when large-format cameras became common, astronomers were once again motivated to generate 3-color composites. CCDs are approximately linear, unlike the roughly logarithmic sensitivity curve of photographic plates, so the natural way to produce an image was to map the logarithm of the pixel intensities to the red, green, and blue brightnesses in the final image (see e.g. Villard and Levay, (2002)). This algorithm produces pictures that look rather like those produced from photographs, but unfortunately it also discards a lot of information.

## 2 Generating Colored Images from CCD Data

Let us consider CCD data taken in three bands which we shall call $`b`$, $`g`$, and $`r`$ (E.g. $`B`$, $`R`$, and $`I`$); we assume that these data have been background-subtracted, if necessary linearised, flatfielded, and scaled appropriately. This scaling is to some extent arbitrary, as different astronomers prefer different color balances, but we have found that a conversion to $`f_\nu `$ (in practice to an AB scale, Oke and Gunn, (1983)) works well.

We desire a mapping to the range $`[0,1]`$ for each of three colors red (R), green (G), and blue (B). The usual algorithm is:

$`R=f(r);G=f(g);B=f(b)`$ (1)
where

$$f(x)=\{\begin{array}{cc}0,\hfill & x<m\text{;}\hfill \\ F(xm)/F(Mm)\hfill & mxM\text{;}\hfill \\ 1\hfill & M<x\text{.}\hfill \end{array}$$

and $`m`$ is the minimum value to display, and $`M`$ the maximum; these are chosen to bring out favoured details. As mentioned in the introduction, a common choice is logarithmic ($`F\mathrm{ln}`$), but linear or square-root stretches are also popular. The actual implementation, at least for integer input data, often uses a lookup table.

This style of mapping simulates one of the unfortunate features of photographic plates, namely that all pixels for which $`r`$, $`g`$, and $`b`$ are $`M`$ appear white in the final image. In a photograph this white color corresponds to a lack of information, but for a CCD it is a choice made by the astronomer. In fact, for *any* non-linear function $`F`$, an object’s color in the composite image depends upon its brightness, even if *none* of $`r`$, $`g`$, or $`b`$ exceed $`M`$.

Fortunately, we can easily avoid these problems. Defining $`I(r+g+b)/3`$ 222It is possible to use a different measure of $`I`$, e.g. a root-mean-square or $`\sqrt{\chi ^2}`$ value., set

$`R`$ $`=rf(I)/I`$
$`G`$ $`=gf(I)/I`$ (2)
$`B`$ $`=bf(I)/I`$
If $`\mathrm{max}_{RGB}\mathrm{max}(R,G,B)>1`$, set $`RR/\mathrm{max}_{RGB}`$, $`GG/\mathrm{max}_{RGB}`$, and $`BB/\mathrm{max}_{RGB}`$; if $`I0`$, set $`R=G=B=0`$. In other words, the intensity is clipped at unity, but the *color* is correct. 333We note that it is often possible to include WCS information in the header of the RBG image; file formats that allow this include JPEG and TIF. Additionally, it is possible to choose a more flexible functional form for $`F`$; following Lupton et al., (1999), we take $`F(x)\text{asinh}(x/\beta )`$, where the softening parameter $`\beta `$ is chosen to bring out desired details — it determines the point at which we change from a linear to a logarithmic transformation ($`\text{asinh}(x)x`$ for $`x1`$; $`\text{asinh}(x)\mathrm{ln}(2x)`$ for $`x1`$) 444A convenient parameterisation is $`f(x)=\text{asinh}(\alpha Q(xm))/Q`$, which allows the user to first set $`Q0`$ and choose the linear stretch $`\alpha `$, and then adjust $`Q`$ to bring out brighter features. In this case, $`Mm+\mathrm{sinh}(Q)/(\alpha Q)`$.. The linear part is used to show faint features in the data, while the logarithmic part enables us to simultaneously illustrate structures such as spiral arms or the cores of star clusters that are significantly brighter. As an additional trick, we find that sometimes a $`\sqrt{\text{asinh}}`$ stretch brings out the details of dust lanes in spiral galaxies; here, as before, the asinh is used to control the dynamic range.

The results of this simple change are dramatic; a given color (i.e. value of $`gb`$ and $`rg`$) is now mapped to a unique color in the RGB image. Fig. 1 shows an SDSS (York et al.,, 2000) g-r-i composite of the galaxies NGC6976, NGC6977, and NGC6978. The reader will notice that the stars, rather than all being white, show a variety of red, yellow, and bluish tints; that the spiral bulges (and the group of ellipticals to the lower right of NGC6977) are a uniform yellowish color, while the spiral disks are blue ($`[OIII]`$ and $`H_\alpha `$ appear in B and R respectively). Furthermore, it is evident that star formation activity in the disks (as indicated by the distinctive bluish color) decreases from left to right.

We may use the well-known color composite of the Hubble Deep Field (HDF; Williams et al, (1996)), Fig. 6, as a pedagogical example. The top row of images were prepared using using Eq. 1, while the bottom row used Eq. 2. The top left panel shows the familiar STScI color image. The bottom left panel use an asinh stretch, and is our preferred representation of the HDF. The top middle panel is our best match to the STScI image and uses a logarithmic stretch; the middle bottom panel uses the same logarithmic function (i.e. identical values of $`m`$ and $`M`$). The righthand pair of panels share a hard linear stretch (it is identical to the linear part of the asinh stretch of the bottom left panel).

Interestingly, the bottom set of panels are pretty similar, with the main difference being that the spiral arms in the lower left corner are rather better defined in the leftmost panel. It is clear that much of the morphological information is carried in the *color* rather than the *intensity*.

The bands ($`f450`$, $`f606`$, and $`f814`$) used in Fig. 6 are not quite the same as those used in the SDSS composites ($`g`$, $`r`$, and $`i`$) but the astrophysical palette is similar, in the sense that the larger HDF galaxies have colors similar to those seen in Fig. 1, with blue disks and yellowish bulges. The fainter objects, however, look unlike *anything* seen in SDSS imaging, with bright red (K-corrected) bulges and, in some cases, bluish disks which must be due to active star formation and the emission of far-uv light.

Fig. 3 shows the globular cluster NGC2419; Fig. 4 shows a color-magnitude diagram of point sources in a 0.6 square degree region surrounding the cluster; the points in this diagram have the *same* color as the corresponding objects in the image; the cluster horizontal- and giant branches stand out.

Finally, this technique is not restricted to optical data; Fig. 7 shows a Chandra X-ray image of the supernova remnant Cas A.

## 3 Detector Saturation

Although CCDs are wonderful detectors, they are not perfect. The SDSS data processing pipeline (Lupton et al.,, 2001) interpolates over saturation trails, which accounts for their absence from these images. Unfortunately, it is not able to correctly recover the counts in the cores of saturated stars, but we *do* know which pixels were involved555If this information isn’t available, all pixels above some threshold may be considered tainted..

Eq. 2 means that every pixel in the image preserves its color, even in the saturated cores of stars. We can avoid undesirable artifacts in the final picture by finding connected pixels which were saturated in at least one band, and replacing them with the average color of the pixels touching the saturated region. Because we do not wish to average over the long bleed-trails produced by bright stars, we only apply this procedure to saturated pixels with intensities $`>M`$. An obvious example of a bright star is seen at the bottom right of Fig. 3; in this case we treated all pixels with at least 10000 counts in any band as ‘saturated’.

## 4 Practical Implementations

The original implementation of these algorithms was included in the SDSS image processing code, but fortunately stand-alone versions are also available. An IDL version is available from http://cosmo.nyu.edu/hogg/visualization/rgb/, and a C version from http://www.sdss.jhu.edu/doc/jpeg.html.

## 5 Conclusions

Color images have traditionally been considered a luxury, but with the techniques espoused in this paper, we have found that they convey an enormous amount of information. For example, star-forming regions have a distinctive color, and it is very clear, looking at images of spiral galaxies, how much the star formation rate varies. We encourage the reader to visit http://www.astro.princeton.edu/$``$rhl/PrettyPictures where, among other objects, all of the NGC galaxies in the SDSS DR1 (Abazajian et al.,, 2003) are presented.

RHL thanks Michael Strauss for helpful comments on an early version of this manuscript.

Funding for the creation and distribution of the SDSS Archive has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Aeronautics and Space Administration, the National Science Foundation, the U.S. Department of Energy, the Japanese Monbukagakusho, and the Max Planck Society. The SDSS Web site is http://www.sdss.org/.

The SDSS is managed by the Astrophysical Research Consortium (ARC) for the Participating Institutions. The Participating Institutions are The University of Chicago, Fermilab, the Institute for Advanced Study, the Japan Participation Group, The Johns Hopkins University, Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, University of Pittsburgh, Princeton University, the United States Naval Observatory, and the University of Washington.

## 6 Figures


text:
# Search for 𝐾⁰_𝑆𝐾⁰_𝑆 in 𝐽/𝜓 and 𝜓(2𝑆) decays

## 1 Introduction

The decay of $`J^{PC}=1^{}`$ charmonium states, like $`J/\psi `$ and $`\psi (2S)`$, to $`K_S^0K_S^0`$ is a CP violating process. However, since CP violations in both the $`K^0\overline{K^0}`$ and $`B^0\overline{B^0}`$ systems have been well established , it is of interest to search for them in other possible processes.

Furthermore, it has been shown  that the $`K_S^0K_S^0`$ system can be used to test the Einstein-Podolsky-Rosen (EPR) paradox versus quantum theory. The space-like, separated coherent quantum system may yield some $`J/\psi K_S^0K_S^0`$ decays if EPR’s locality is correct while quantum theory forbids this process.

MARK-III searched for decays of $`J/\psi K_S^0K_S^0`$ with 2.7 million $`J/\psi `$ events . No signal was observed, and the upper limit on the decay branching ratio was determined to be $`(J/\psi K_S^0K_S^0)<5.2\times 10^6`$ at the 90% confidence level (C. L.). The same search for $`\psi (2S)`$ decays has not been performed before.

In this Letter, we report on a search for $`J/\psi K_S^0K_S^0`$ using a sample with 20 times more statistics than before, and on the first search for $`\psi (2S)K_S^0K_S^0`$. The data samples used for the analyses are taken with the Beijing Spectrometer (BESII) detector at the Beijing Spectrometer (BEPC) storage ring at a center-of-mass energies corresponding to $`M_{J/\psi }`$ and $`M_{\psi (2S)}`$. The data samples contain $`(57.7\pm 2.7)\times 10^6`$ $`J/\psi `$ events and $`(14\pm 0.7)\times 10^6`$ $`\psi (2S)`$ events, as determined from inclusive hadronic decays .

BES is a conventional solenoidal magnet detector that is described in detail in Ref. ; BESII is the upgraded version of the detector . A 12-layer vertex chamber (VC) surrounding the beam pipe provides trigger information. A forty-layer main drift chamber (MDC), located radially outside the VC, provides trajectory and energy loss ($`dE/dx`$) information for charged tracks over $`85\%`$ of the total solid angle. The momentum resolution is $`\sigma _p/p=0.017\sqrt{1+p^2}`$ ($`p`$ in $`\text{GeV}/c`$), and the $`dE/dx`$ resolution for hadron tracks is $`8\%`$. An array of 48 scintillation counters surrounding the MDC measures the time-of-flight (TOF) of charged tracks with a resolution of $`200`$ ps for hadrons. Radially outside the TOF system is a 12 radiation length, lead-gas barrel shower counter (BSC). This measures the energies of electrons and photons over $`80\%`$ of the total solid angle with an energy resolution of $`\sigma _E/E=22\%/\sqrt{E}`$ ($`E`$ in GeV). Outside of the solenoidal coil, which provides a 0.4 Tesla magnetic field over the tracking volume, is an iron flux return that is instrumented with three double layers of counters that identify muons of momentum greater than 0.5 GeV/$`c`$.

## 2 Monte Carlo

A Monte Carlo (MC) simulation is used for the determination of mass resolutions and detection efficiencies. For the signal channels, $`J/\psi K_S^0K_S^0`$ and $`\psi (2S)K_S^0K_S^0`$, the angular distribution of one $`K_S^0`$ is generated with the same distribution as for the CP allowed channels $`J/\psi K_S^0K_L^0`$ and $`\psi (2S)K_S^0K_L^0`$, namely as $`\mathrm{sin}^2\theta `$, where $`\theta `$ is the polar angle of the $`K_S^0`$ in the laboratory system. Only $`K_S^0\pi ^+\pi ^{}`$ is generated, and the $`K_S^0`$ is allowed to decay in the detector according to its lifetime. The simulation of the detector response is done using a Geant3 based program with detailed consideration of the detector performance (such as dead electronic channels). Reasonable agreement between data and Monte Carlo simulation has been observed in various channels tested, including $`e^+e^{}(\gamma )e^+e^{}`$, $`e^+e^{}(\gamma )\mu ^+\mu ^{}`$, $`J/\psi p\overline{p}`$ and $`\psi (2S)\pi ^+\pi ^{}J/\psi ,J/\psi \mathrm{}^+\mathrm{}^{}`$ $`(\mathrm{}=e,\mu )`$. For this study, 10,000 events for each signal channel are produced.

## 3 Event selection

For the decay channels of interest, the candidate events are required to satisfy the following selection criteria:

1. The number of charged tracks is required to be four with net charge zero. Each track should satisfy $`|\mathrm{cos}\theta |<0.8`$, where $`\theta `$ is the polar angle in the MDC, and should have a good helix fit so that the error matrix of the track fitting is available for secondary vertex finding.
2. The tracks are assumed to be either $`\pi ^+`$ or $`\pi ^{}`$. The higher momentum positive track and the lower momentum negative track are assumed to come from one of the $`K_S^0`$ decays, and the remaining two tracks are assumed to come from the other $`K_S^0`$ decay. The intersection of the track-pairs near the interaction point are determined and are regarded as the $`K_S^0`$ vertices.
3. Each $`K_S^0`$ candidate is required to have a decay length in the transverse plane ($`L_{xy}`$) greater than 3 mm and an invariant mass within twice the mass resolution as estimated from Monte Carlo simulations (the mass resolution in $`J/\psi `$ decays is 7.1 MeV/$`c^2`$ and in $`\psi (2S)`$ decays 7.9 MeV/$`c^2`$).

For $`\psi (2S)`$ decays, backgrounds from $`\psi (2S)\pi ^+\pi ^{}J/\psi `$ with $`J/\psi `$ decaying into two charged tracks are removed by further requiring the mass recoiling against the lower momentum positive and negative tracks, $`m_{rec}`$, to be less than 3.0 GeV/$`c^2`$. The comparison of the recoil mass distributions of signal and background events from Monte Carlo simulation is shown in Figure 1. It can be seen that this requirement is very efficient in removing background, and the loss of efficiency is small ($``$ 4%).

Another possible background is from $`\psi (2S)\gamma \chi _{cJ}`$ $`(J=0,2)`$ with $`\chi _{cJ}K_S^0K_S^0`$. Figure 2 shows $`K_S^0K_S^0`$ missing momentum distributions for the signal and background channels, all simulated by Monte Carlo. We make an additional requirement $`p_{miss}<0.1`$ GeV/$`c`$ for $`\psi (2S)`$ selection. This removes 94% of $`\chi _{c2}`$ background and almost all the $`\chi _{c0}`$ background, while the efficiency for $`\psi (2S)K_S^0K_S^0`$ is about 92% according to Monte Carlo simulation. For $`J/\psi `$ decays, since $`\eta _cK_S^0K_S^0`$ is forbidden by parity conservation, no requirement is applied.

The momenta of the two $`K_S^0`$ tracks for events passing the above selection are plotted in Figures 3 and 4, for $`J/\psi `$ and $`\psi (2S)`$ decays, respectively, together with the regions predicted for signal events by the Monte Carlo simulation. The two circles correspond to

$$\sqrt{(p_{K_S^0}^ap_0)^2+(p_{K_S^0}^bp_0)^2}<n\sigma ,$$

where $`p_{K_S^0}^a`$ and $`p_{K_S^0}^b`$ are the momenta of the two $`K_S^0`$ candidates in the event; $`\sigma `$ is the $`K_S^0`$ momentum resolution, which is 28 MeV/$`c`$ in $`J/\psi `$ decays and 34 MeV/$`c`$ in $`\psi (2S)`$ decays, as determined from Monte Carlo simulations; and $`p_0`$ is the expected $`K_S^0`$ momentum in $`J/\psi `$ (1.466 GeV) or $`\psi (2S)`$ (1.775 GeV) two-body decays. The inner circles correspond to $`n=1`$, and the outer circles are for $`n=2`$. It can be seen that there are no candidate events in either $`J/\psi `$ or $`\psi (2S)`$ decays within the one $`\sigma `$ circle, and there is only one candidate in each case lying between the one and two $`\sigma `$ circles. Using the 2$`\sigma `$ circle for our event selection, one event is obtained each for $`J/\psi K_S^0K_S^0`$ and for $`\psi (2S)K_S^0K_S^0`$. The corresponding efficiencies, $`\epsilon _{MC}`$, are 20.74% and 19.18% for $`J/\psi `$ and $`\psi (2S)`$ decays, respectively, as estimated from Monte Carlo samples.

The backgrounds remaining in both $`J/\psi `$ and $`\psi (2S)`$ decays after the $`K_S^0`$ requirement are studied with Monte Carlo simulations. The main background channel in $`J/\psi `$ decays is from $`J/\psi \overline{K}^{}(892)^0K^0+c.c.`$ is estimated to be $`1.0\pm 0.4`$ events, which is obtained using a Monte Carlo sample six times as large as the $`J/\psi `$ sample. The main backgrounds channels in $`\psi (2S)`$ decays are $`\psi (2S)\gamma \chi _{c2}`$, $`\chi _{c2}K_S^0K_S^0`$ and $`\psi (2S)\overline{K}^{}(892)^0K^0+c.c.`$. A Monte Carlo simulation estimates the total background to be $`0.67\pm 0.25`$ events, with about 80% due to $`\chi _{c2}K_S^0K_S^0`$. The errors on the estimated backgrounds are from the uncertainties of the branching ratios used and the limited statistics of the Monte Carlo samples.

The two candidate events are investigated further, and it is found that the $`J/\psi K_S^0K_S^0`$ candidate is most likely a $`J/\psi \overline{K}^{}(892)^0K^0+c.c.`$, with $`K^0K^{}\pi ^+`$ misidentified as $`K_S^0\pi ^+\pi ^{}`$ because no particle identification is used in the event selection, while the $`\psi (2S)K_S^0K_S^0`$ candidate looks like a $`\psi (2S)\gamma \chi _{c2}`$, $`\chi _{c2}K_S^0K_S^0`$ event with an isolated photon with BSC energy of 126 MeV.

## 4 Trigger efficiency

Due to the long decay length of the high momentum $`K_S^0`$ and the requirement of hits in the VC, the trigger efficiency of $`K_S^0K_S^0`$ events is lower than for normal hadronic events. Since the trigger system is not included in the Monte Carlo simulation, the trigger efficiency is measured with data by comparing events within and beyond the outer radius of the VC using $`J/\psi (\psi (2S))K_S^0K_L^0`$ events, which yields trigger efficiencies of $`(86.7\pm 0.9)\%`$ and $`(81.3\pm 1.9)\%`$ for a single $`K_S^0`$ in $`J/\psi `$ and $`\psi (2S)`$ decays, respectively. For $`K_S^0K_S^0`$ events, the trigger efficiency, $`\epsilon _{trig}`$, is calculated to be $`(98.2\pm 0.2)`$% in $`J/\psi `$ decays and $`(96.5\pm 0.7)`$% in $`\psi (2S)`$ decays.

## 5 Secondary Vertex Finding

The secondary vertex finding algorithm is checked using $`J/\psi K^{}(892)\overline{K}+c.c.`$ events. It is found that the Monte Carlo simulates data fairly well, and extrapolating the difference between data and Monte Carlo simulation to the $`K_S^0`$ momentum range under study and correcting by the polar angle dependence of the efficiency, correction factors of $`(96.4\pm 3.1)\%`$ and $`(98.1\pm 4.0)\%`$ are obtained for a single $`K_S^0`$ in $`J/\psi `$ and $`\psi (2S)`$ decays, respectively. For $`K_S^0K_S^0`$ decays, the correction factors and errors to the Monte Carlo efficiencies, $`\epsilon _{2nd}`$, for $`J/\psi K_S^0K_S^0`$ and $`\psi (2S)K_S^0K_S^0`$ are $`(92.9\pm 4.5)\%`$ and $`(96.2\pm 5.8)\%`$, respectively.

## 6 Systematic error

The total systematic error on the branching ratio measurement comes from all sources listed in Table. 1. The simulation of the tracking efficiency agrees with data within 1-2% for each charged track as measured using channels like $`J/\psi \mathrm{\Lambda }\overline{\mathrm{\Lambda }}`$ and $`\psi (2S)\pi ^+\pi ^{}J/\psi `$, $`J/\psi \mu ^+\mu ^{}`$, and 8% is taken conservatively as the systematic error for the channel of interest.

The Monte Carlo simulated $`K_S^0`$ mass and momentum resolutions agree with those of data within statistical uncertainties, as has been checked with $`J/\psi `$ and $`\psi (2S)K_S^0K_L^0`$ samples . The requirement that the $`K_S^0`$ mass and momentum be within two standard deviations introduces systematic uncertainties at the 5-6% level, dominated by the statistical precisions of the comparisons between data and Monte Carlo simulation.

The systematic errors on the $`p_{miss}`$ and the recoil mass requirements for the $`\psi (2S)`$ sample depends on the simulation of the momentum of the charged tracks and are already included either in the tracking or in the systematic error quoted for the $`K_S^0`$ mass and momentum resolutions. They are not further considered here.

The systematic errors on the total numbers of $`J/\psi `$ and $`\psi (2S)`$ events are taken as 4.7% and 5.0%, respectively, and are measured using inclusive hadronic events with four charged hadrons in the final state  for $`J/\psi `$ and using inclusive hadrons for $`\psi (2S)`$ . The systematic error on the branching ratio used, $`(K_S^0\pi ^+\pi ^{})`$ is obtained from the Particle Data Group  directly. Adding all the systematic errors in quadrature, the total systematic errors are 12.0% and 12.7% for $`J/\psi `$ and $`\psi (2S)`$ decays, respectively.

## 7 Results

Since the observed numbers of events agree with the expected background levels for the channels under study, upper limits are conservatively set by not subtracting background and taking the events in MC predicted region as signal. With one observed event, the upper limit on the number of events is 4.74 at the 95% confidence level.

The upper limits on the branching ratios of $`J/\psi `$ and $`\psi (2S)K_S^0K_S^0`$ are calculated with

$$(RK_S^0K_S^0)<\frac{n_{UL}^{obs}/(\epsilon _{MC}\epsilon _{trig}\epsilon _{2nd})}{N_R(K_S^0\pi ^+\pi ^{})^2},$$

where $`n_{UL}^{obs}`$ is the upper limit of the number of observed events, and $`N_R`$ is the total number of resonance $`R`$ events. Using the numbers listed in Table 2, one obtains the upper limits on the branching ratios at the 95% C. L.:

$$(J/\psi K_S^0K_S^0)<1.0\times 10^6,$$

$$(\psi (2S)K_S^0K_S^0)<4.6\times 10^6,$$

where the systematic errors are included by lowering the efficiencies by one standard deviation.

## 8 Summary

The CP violating processes $`J/\psi K_S^0K_S^0`$ and $`\psi (2S)K_S^0K_S^0`$ are searched for using the BESII 58 million $`J/\psi `$ event and 14 million $`\psi (2S)`$ event samples. The upper limits on the branching ratios are determined to be $`(J/\psi K_S^0K_S^0)<1.0\times 10^6`$ and $`(\psi (2S)K_S^0K_S^0)<4.6\times 10^6`$ at the 95% C. L. The former is much more stringent than the previous MARK-III measurement , and the latter is the first search for this channel in $`\psi (2S)`$ decays. Current bounds of the production rates still beyond the sensitivity for testing the EPR paradox .

## Acknowledgments

We wish to thank Dr. H. B. Li for helpful comments. The BES collaboration thanks the staff of BEPC for their hard efforts. This work is supported in part by the National Natural Science Foundation of China under contracts Nos. 19991480, 10225524, 10225525, the Chinese Academy of Sciences under contract No. KJ 95T-03, the 100 Talents Program of CAS under Contract Nos. U-11, U-24, U-25, and the Knowledge Innovation Project of CAS under Contract Nos. U-602, U-34(IHEP); by the National Natural Science Foundation of China under Contract No. 10175060 (USTC); and by the Department of Energy under Contract No. DE-FG03-94ER40833 (U Hawaii).


text:
# An Unbinned Goodness-of-Fit Test Based on the Random Walk1footnote 11footnote 1University of Cincinnati preprint # UCHEP-03-02

## I Introduction

This search for an unbinned goodness-of-fit test has been motivated by the widespread use of unbinned maximum likelihood fitting for determining $`CP`$-violating parameters at Belle. While there are many cross-checks to insure that there are no spurious signals and biases, the fits tend to be complicated and not very transparent. They often involve probability density functions (PDF’s) that differ with every event, based on measured quantities that add dimensions to the data that are not explicit in the fits. As there is no widely accepted unbinned goodness-of-fit test that applies to such fits, testing for statistical consistency of results has been uneven. The tests that have been done, resorting to binned $`\chi ^2`$ or toy Monte Carlo, have their place but have not been entirely satisfactory in addressing the question.

A common technique of unbinned tests involves first transforming the measured quantities to a variable in which the null hypothesis has a uniform distribution, where the PDF is flat, and then to test this “flattened” distribution for consistency with uniformity. There exists a variety of tests for uniformity, but most are not readily extended to multidimensional data, and they do not address compound hypotheses. A review of methods is given in aslandurham .

In this report, we explore a test statistic that is based on the two-dimensional Random Walk. To begin, its distribution in the case of a flat PDF is discussed. The ensemble distribution is then found for several alternate hypotheses, and the rejection power is calculated for comparison with other goodness-of-fit tests. As the aim of a goodness-of-fit test as it would be applied at Belle is to test the validity of the parametrization used in fitting, it is also important to examine how the test is modified under compound hypotheses. The discussion is thus expanded to include data which are fitted to determine one or more parameters. Finally, we discuss the possibility of extending to multidimensional data.

## II Random Walk as a Test of Flatness

A dataset consisting of $`N`$ measurements of the one-dimensional quantity $`x`$ lying in the interval $`[0,1]`$ may be mapped trivially to points on a unit circle with polar angle $`\varphi `$ on the interval $`[0,2\pi ]`$, so that each point is considered to be a unit vector with direction defined by $`\varphi `$. If the PDF in $`x`$ is flat, the vector sum of the corresponding unit vectors in two dimensions corresponds to the net displacement, $`D`$, after a two-dimensional Random Walk of $`N`$ steps with unit step size. For sufficiently large $`N`$, this distribution converges to a well-known form (Rayleigh, 1888) and the distribution in $`D^2`$ is an exponential decay with mean equal to $`N`$ . We take $`D^2/N`$ as the test statistic. A deviation of the root distribution from the hypothesis will result in a bias of the ensemble distribution of this test statistic away from the origin. This statistic is mathematically equivalent to the first order term in the Fourier series that describes the distribution of the data:

$`(k=1)`$ $`=`$ $`{\displaystyle _0^{2\pi }}𝑑\varphi {\displaystyle \underset{j=1}{\overset{N}{}}}e^{ik\varphi }\delta (\varphi \varphi _j)`$
$`=`$ $`{\displaystyle \underset{j=1}{\overset{N}{}}}e^{i\varphi _j}`$
where one can see that $`D^2|(1)|^2`$. One would expect this distribution to be most sensitive to an overall imbalance of the PDF in generally opposite $`\varphi `$ directions. To obtain sensitivity to higher order differences, one could thus take successively higher order terms in the series, for $`k=2,\mathrm{}`$. In practice it may not be useful to examine terms above $`k=3`$. In this study we look at $`k=1`$ ($`d=1`$) and define $`K_k\frac{|(k)|^2}{N}`$. What we have defined as $`K_1`$ appears in the review of D’Agostino and Stephensstatistic  as $`R`$ in the context of the Von Mises test, a test for uniformity on a circle.

## III Flat PDF

As mentioned above, the $`K_1`$ distribution for a flat PDF converges rapidly to an exponential with a decay constant of unity. Figure 1 (top row) shows the distributions in $`K_1`$ for ensembles of randomly generated experiments containing $`N=10`$, 100, and 1000 events. Each of the three distributions is fitted via binned maximum likelihood to an exponential form. The fitted inverse decay constants (“slopes”) are $`0.992\pm 0.010`$, $`1.008\pm 0.033`$, and $`1.039\pm 0.049`$, respectively, in excellent agreement with the expectation.

To evaluate rejection power, these distributions may be compared with those obtained for PDF’s that are not flat. The alternative hypotheses used in a study by Aslan and Zech Zech  provide a convenient range of function types and allow for a direct comparson with the range of tests reviewed in their work. In that paper the rejection power of the alternative hypothesis was defined as one minus the probability for an error of the second kind, given a criterion that yields a 5% significance for the null hypothesis. Since in this case the null hypothesis gives an exponential distribution with unit decay constant, the 5% criterion is $`K_1>3.0`$. Ensembles of experiments were generated for each of three functions used in Ref. Zech :

$`𝒜_1(X)`$ $`=`$ $`0.3+1.4X`$ (2)
$`𝒜_2(X)`$ $`=`$ $`0.7+0.3[n_2e^{64(X0.5)^2}]`$ (3)
$`𝒜_3(X)`$ $`=`$ $`0.8+0.2[n_3e^{256(X0.5)^2}]`$ (4)
where the $`n_i`$ are normalization constants for the associated Gaussians. All functions are defined in the interval $`[0,1]`$. The resulting $`K_1`$ distributions for $`𝒜_1`$ are shown in Figure 1 (bottom row). The values for rejection power are summarized in Table 1. For comparison, the values for the $`\chi ^2`$ method ($`N=100`$) given by Ref. Zech  are approximately 0.81, 0.85, and 0.81, respectively, so our method is comparable in power, at least in the case of these three functions.

In order to apply this method as a goodness-of-fit test for non-uniform null hypotheses the PDF, $`f(X)`$, must first be transformed to a “flat” variable, $`Y`$, where the probability distribution is flat. To form a uniform null hypothesis on a circle one could, for example, construct $`Y`$ as :

$`Y_i`$ $`=`$ $`2\pi {\displaystyle _X_{}^{X_i}}f(X)𝑑X`$ (5)
where the integer subscript $`i`$ denotes the $`i^{th}`$ data point and $`X_{}`$ is the lowest possible value of $`X`$.

## IV Compound hypotheses

The examples considered thus far have been ones where no parameter fitting has occurred. While this has been an instructive exercise, it has limited application, as most measurements in particle physics involve the fitting of measured distributions to determine shapes and to derive some physics quantity or conclusion. We now look at compound hypotheses.

In evaluating rejection of alternative hypotheses via toy MC in the compound case, it is important that the fitting process be integrated into the evaluation procedure. Consider a data set $`\{\varphi _i\}`$ where the PDF is assumed to be parametrizable as $`f(\varphi ;\alpha )`$ and the unbinned likelihood is maximum for $`\alpha =\alpha _{max}`$. The data are then flattened assuming the PDF is $`f(\varphi ;\alpha _{max})`$, and the associated $`K_1`$ is evaluated. The confidence level of this $`K_1`$ value may then be found by referencing the ensemble distribution of $`K_1`$ when the true PDF is $`f(\varphi ;\alpha _{max})`$, and each experiment of the ensemble is treated as data, fitted and flattened according to the fit.

This procedure was used to evaluate rejection power for pairs of similarly shaped PDF’s. Here we show one such result, for the hypothesis $`n_4(\alpha )e^{10X\alpha }`$, where $`n_4`$ is a normalization constant, the measured quantity is $`X`$, and experiments are fitted for $`\alpha `$. The alternative PDF was the linear form $`f(X)=2(1X)`$. Experiments were generated according to the alternative PDF (A), and each was fitted to the hypothesis. The mean maximum likelihood value of $`\alpha `$ was approximately 4.7. Ensembles (B) were generated according to the hypothesis, with $`\alpha =4.7`$, and fitted in the same way. The 5% confidence criterion on $`K_1`$ for (B) and acceptance of this criterion for (A) were estimated by counting (Figure 2). The rejection powers were found to be 28% and 99% for $`N=100`$ and $`N=1000`$, respectively. For comparison we also calculated by the same procedure the rejection of the $`\chi ^2`$ test, using 20 bins in the interval  and found powers of 13% and 100%, respectively.

We also examined the two-dimensional distribution of fitted $`\alpha _{max}`$ values $`vs.`$ $`K_1`$. Any dependence of the test on the fitted rather than underlying parameter value reduces its utility as a goodness-of-fit test; for example, the maximum likelihood value, $`_{max}`$, is not usable as a goodness-of-fit statistic because it depends strongly on the fitted parameter value(s) $`\alpha _{max}`$ – for a certain class of fitting functions, the correlation is 100%likelihood . Figure 3 shows scatter plots of $`\alpha _{max}`$ and $`K_1`$, where the data were generated with $`N=100`$ and the generated and fitted forms are those from the example of Figure 2. There appears to be no strong dependence.

In any determination of rejection power with a compound hypothesis, it is necessary to determine the distribution of $`K_1`$ for the correct hypothesis. It does not appear that there is a simple ansatz as in the case of binned least squares fitting, where the chisquare converges to a chisquare distribution with the number of degrees of freedom reduced by one unit for each linear fitted parameter. We study this question empirically by generating MC ensembles for a variety of shapes. Each ensemble was generated according to the fitted functional form with parameter value(s) fixed. Each experiment was fitted with parameter(s) floating, and the $`K_1`$ value was obtained from the data flattened according to the best fit. The distribution of resultant $`K_1`$ values for each ensemble was fitted for the decay constant, assuming an exponentially decaying form. Ensembles with $`N=10`$, $`N=100`$, and $`N=1000`$ were generated. The results are summarized in Table 2. There are several notable features. First, while all of the $`K_1`$ distributions had a decaying form, as one might expect, and a fit that converged, not all yielded good fits; the exponential form is not preserved under compound hypotheses. Secondly, all inverse decay constants are greater than unity, indicating that the $`K_1`$ distribution moves toward zero with fitting. This is not suprising; fitting identifies for each experiment the shape that is “closest” to the data, giving in general a better goodness-of-fit than the generator shape. Finally, there is no obvious pattern in the value of the decay constant with number of floated parameters. However, it is seen that for a given PDF and set of fitted parameters, the shape of the $`K_1`$ distribution shows remarkably little change as $`N`$ is changed by two orders of magnitude.

## V Extension to multidimensional data: speculation

Our goal in this investigation has been to arrive at a multidimensional unbinned goodness-of-fit test, one that has rejection power in all dimensions, not just in one-dimensional projections, for multidimensional data. Many unbinned tests depend on the integrated sum of or spacings between neighboring data points, quantities which are not well-defined when extended to more than one dimension. Although the $`K_1`$ statistic does not have this property, it is yet to be determined whether there exists an extension that is fully multidimensional; for example, in two dimensions, two components each mapped to a circle corresponds to a data space that is the surface of a toroid, for which there is no obvious nontrivial vector sum that maps to the Random Walk. A fully general extension to multidimensional data will additionally require a flattening algorithm and provisions for data spaces of arbitrary shape. We will continue to explore the possibilities for extending $`K_1`$ for use with multidimensional data.

## VI Summary

We have explored an unbinned goodness-of-fit test for data in one dimension that is based on the mapping of flattened distributions to a two-dimensional random walk. This method is truly binning-free and scale-independent, and the ensemble distribution for the null hypothesis is well-defined. For a compound hypothesis we specify a procedure to determine the ensemble distribution of the test statistic via Monte Carlo so that rejection power may be readily determined. The distribution is found for several different parametrized forms and shown to be largely independent of statistics. We examine several samples for dependence between the test statistic and fitted parameter values, and find no evidence of any. The rejection power for alternate hypotheses is demonstrated for a few examples and is found to be comparable to that of the chisquare method.

###### Acknowledgements.

The author would like to thank R. Cousins, G. Zech, and B. Yabsley for useful discussions and suggestions, and the organizers of PHYSTAT 2003 for a stimulating and interesting conference. This work is supported by Department of Energy grant #DE-FG02-84ER40153.


text:
# 1 Introduction

## 1 Introduction

As the experimental data on the large-scale structure and early history of the universe improve, there has been increased interest in seeing how ideas from string theory may fit in to understanding cosmology at the earliest of times. On the whole, efforts in this direction have focused on either finding a home for inflation in string theory (for instance, ) or replacing it with something quite different but string-inspired (for instance, ). Important as it is to establish in a fundamental theory of gravity either a sound underpinning for inflation or a solid alternative to it, my focus here will be instead on the possibility of string physics playing an important role just after inflation has ended. The characteristic scale of string theory is usually closely to the Planck mass, considerably higher than a typical inflaton mass scale. However, it turns out that objects known as “tensionless strings”  can be produced copiously from coherent motion of the inflaton field. Without significant fine-tuning, it is possible that creation of such strings can extract a significant fraction of energy from the inflaton, thus providing an alternative to preheating .

String creation is interesting in its own right, but on the whole it is ill-understood. The best studied examples of it are decaying D-brane solutions , where boundary CFT methods are available. We will rely instead on estimates from effective field theory, employing a method of steepest descent to extract Bogliubov coefficients from approximate wave equations for massive string modes. This method was developed in , and its application to preheating was considered in . It would be highly desirable to improve our ability to compute string creation amplitudes: this is an important feature of time-dependent backgrounds in string theory, and it is relevant to almost all the stringy proposals for the physics of the early universe.

The organization of this paper is as follows. In section 2 we will briefly summarize the results of the steepest descent analysis of string creation. In section 3 we will explain in some detail how this analysis may be applied to provide an alternative to preheating. In section 4 we will remark on a special background in which a worldsheet treatment of closed string creation may be possible.

The results summarized here are based mainly on , except that section 4 is based on ongoing work with J. Friess and I. Mitra.

## 2 Estimating string creation

The density of states in tree-level string theory is exponential:

(1)
where $`T_H`$ is the Hagedorn temperature and $`\gamma `$ depends on details of the string theory in question. The quantity $`c_{}`$ is the central charge for transverse modes: $`c_{}=12`$ for the type II superstring. The exponential density of states (1) is in contrast with the power-law behavior one encounters for point-particle theories. The states that dominate (1) are highly excited strings. String interactions are not taken into account in (1). $`T_H`$ appears to be a temperature at which some phase transition occurs. Above the phase transition strings themselves may not be a good description of the physics. There is an extensive literature on this so-called Hagedorn transition (see for example ), and string interactions do enter into the discussion in important ways. For the purposes of the present discussion, we will not need to consider interactions explicitly.

The on-shell constraints for quantum states of a string include $`L_0|\psi _{\mathrm{phys}}=0`$, which leads to an equation roughly of the form

(2)
where $`\omega (t)^2`$ depends on the excitation level and the momentum: in a flat background, $`\omega ^2=k^2+N/\alpha ^{}`$, where $`N`$ is the level (up to factors of $`2`$ that depend somewhat on details). In a cosmological background, a reasonable first guess is

(3)
where $`a(t)`$ is the scale factor and $`t`$ is conformal time. One would of course like to improve upon (2): this is one place where input from a worldsheet treatment is needed. The form (3) is predicated on working with the string frame metric, where the parameter $`\alpha ^{}`$ is a constant length squared.

The steepest descent method, developed in this context in , allows us to crudely estimate the Bogliubov coefficients for (2): if $`|\beta |1`$, then

(4)
where $`t_{}=ri\mu `$ is the location of the zero of $`\omega (t)^2`$ in the lower half plane which is closest to the real axis, and $`r`$ and $`\mu `$ are real. The last of the three factors in (4) is the most important: it leads to exponential suppression of $`|\beta |^2`$ when $`\omega `$ is large. Assuming that $`\omega (t)`$ is well-approximated by an elliptical arc over the contour from $`r`$ to $`t_{}`$, one has

(5)
From this one can estimate the total number of strings created:

| $`N_{\mathrm{tot}}`$ | $`={\displaystyle ^{\mathrm{}}}𝑑\omega {\displaystyle \frac{dN}{d\omega }}|\beta |^2{\displaystyle ^{\mathrm{}}}𝑑\omega e^{\omega /T_H}e^{\pi \mu \omega }.`$ |
| --- | --- |

(6)
Evidently, this total number converges precisely if

(7)
If this bound is violated, it’s somewhat like passing to a temperature above the Hagedorn temperature: clearly, interactions and/or back-reaction must then be involved in rendering the total energy passing into strings finite.

The punch-line of the next section will be that, for non-perturbative strings, it is possible without appreciable fine-tuning of parameters to violate (7) just after inflation. If this happens, then strings must indeed play an important role in reheating. Moreover, it would be highly excited strings that would make the dominant contribution, so the effects of strings would be quite unlike anything in a point-particle theory—hence characteristic of string theory.

Needless to say, it would be very desirable to learn what happens in time-dependent backgrounds when the number or energy of strings created diverges. Decaying D-branes may be the simplest laboratory in which to attack this question: see for example .

## 3 Extracting energy from inflaton oscillations

Let us start with a brief overview of the theory of preheating , which will be a useful point of comparison when considering post-inflationary string creation.

* At the end of inflation, the universe is cold and empty.
* The inflaton $`\varphi `$ is oscillating around its minimum with an amplitude assumed to be comparable to the Planck scale.
* These coherent oscillations may die off via perturbative production of particles. This is ordinary reheating.
* Preheating can occur if another boson has a coupling $`g\varphi ^2\chi ^2`$ to the inflaton.
* The effective mass of $`\chi `$ varies: $`m_{\mathrm{eff}}^2=m_0^2+g\varphi ^2`$. This means parametric resonance can occur and give exponentially growing occupation numbers for $`\chi `$.
* Preheating is “efficient” iff $`m_{\mathrm{eff}}^2`$ varies by a large factor from its minimum to its maximum: $`m_0^2`$ needs to be a small contribution to the mass.

The last point appears contrived, especially when it is stated more precisely, as we will do shortly. “Efficient” in this context means that a significant fraction of the energy in the coherent inflaton oscillations should be extracted via the parametric resonance. We may then ask, if string states could be produced in a similar way, could they play a role in reheating?

We have already noted in the introduction that the string scale is usually only slightly smaller than the four-dimensional Planck scale, $`M_{\mathrm{Pl},4}=(8\pi G_N)^{1/2}10^{18}\mathrm{GeV}`$, while a typical value for the inflaton mass (in power-law inflation) is $`m_{\mathrm{inflaton}}10^{13}\mathrm{GeV}`$. It will not come as a surprise that this prevents copious production of fundamental strings, unless we make somewhat unnatural assumptions. Although our eventual aim is to consider production of non-perturbatively constructed strings which are tensionless at some point in moduli space, it will be instructive first to consider the production of fundamental strings in slightly more detail; also, this discussion will tie in with the special background discussed in section 4.

The first goal in assessing a model of preheating is to find the regime of parameters where preheating is efficient. The analogous goal here is to find where the bound (7) is first violated. When this happens, the particle occupation numbers of the individual string modes should still be small: if they are not, then according to (6), the number of strings created will diverge exponentially. Thus we have no need of parametric resonance and can continue to use the steepest descent method.

The onshell condition in string metric is, approximately, $`\omega ^2=k^2+m^2|g_{tt}^{str}|`$, where $`t`$ is conformal time. The string metric is related to the Einstein metric by

(8)
where $`\phi `$ is the dilaton and $`\gamma `$ is a constant of order unity. Thus

(9)
The dilaton is probably not the same field as the inflaton; however, it seems reasonable to assume some finite overlap, so when the inflaton is oscillating at the end of inflation, the dilaton is too. Thus $`\omega (t)^2`$ oscillates between slowly varying extremes:

| $`\omega (t)^2`$ | $`={\displaystyle \frac{\omega _{\mathrm{max}}^2+\omega _{\mathrm{min}}^2}{2}}+{\displaystyle \frac{\omega _{\mathrm{max}}^2\omega _{\mathrm{min}}^2}{2}}\mathrm{cos}\mathrm{\Omega }t`$ |
| --- | --- |
| $`\omega _{\mathrm{min}}^2`$ | $`=k^2a(t)^2+e^{\gamma \phi _{\mathrm{min}}/M_{\mathrm{Pl},4}}{\displaystyle \frac{N}{\alpha ^{}}}`$ |
| $`\omega _{\mathrm{max}}^2`$ | $`=k^2a(t)^2+e^{\gamma \phi _{\mathrm{max}}/M_{\mathrm{Pl},4}}{\displaystyle \frac{N}{\alpha ^{}}}.`$ |

(10)
The zeroes of $`\omega (t)^2`$ are

(11)
and it is plausible to assume that $`t_{}=t_{,0}^{}=ri\mu `$ makes the biggest contribution to string production. String production is finite if

(12)
which translates approximately to

(13)
where $`\tau _{\mathrm{min}}`$ is the minimum string tension as measured in Einstein frame.

If the bound (12) is violated, then what must happen is that a significant fraction of the energy of the inflaton oscillations passes into excited string states, and the back-reaction slows down the inflaton. This would correspond to the idea of efficient preheating. Note however that $`\omega _{\mathrm{min}}\sqrt{\tau _{\mathrm{min}}}`$ and $`\mathrm{\Omega }m_{\mathrm{inflaton}}`$, so for (13) to occur, the string tension as measured in Einstein frame must dip temporarily below $`m_{\mathrm{inflaton}}`$. This is possible, but it seems like asking a lot in conventional string compactifications. The upshot, then, is that dimensional analysis didn’t fail, and we have proposed a contrived scenario.

Let’s compare with conventional preheating, where it is useful to define “Mathieu parameters” $`A`$ and $`q`$ as follows:

(14)
A “broad resonance” occurs in the standard theory when $`2q<A<2q+\sqrt{bq}`$ for some order one constant $`b`$: this leads to efficient particle production. On the other hand, efficient string production occurs when $`2q<A<2xq`$ for some $`x>1`$. So string production occurs in a qualitatively bigger wedge of parameter space, albeit still a small one.

Another problem with proposing that fundamental strings might have been created after inflation is that they can wind around the extra dimensions. In phenomenologically interesting string compactifications, there usually is a non-trivial first homotopy group on the compactification manifold, so such windings can be topologically stable. Moreover, the density of states of wound strings is smaller only by roughly a constant factor compared to unwound strings, and this constant factor is at most a couple of orders of magnitude in conventional setups. Wound strings would act as dark matter, or perhaps they could even have visible sector effects like fractional charge. Just the dark matter consideration is discouraging: a naive estimate suggests that wound strings would overclose the universe today by something close to a factor of $`10^{17}`$.

Considering non-perturbatively constructed strings in place of fundamental strings alleviates all the difficulties so far discussed. These strings would arise as branes wrapped on a shrinking cycle in the extra dimensions. Consider for example D3-branes on a shrinking $`S^2`$. The $`S^2`$ may become small only at a point on the compactification manifold. This avoids the difficulty of having wound strings, because there is a strong energy barrier against wrapped branes wandering away from the location of the shrinking cycle. The parameter controlling the size of the $`S^2`$ is in fact complex: it is

(15)
where $`J_2`$ is the Kahler form on the extra dimensions (assumed to form a complex manifold) and $`B_2`$ is the Neveu-Schwarz two-form potential of string theory. Starting from the DBI action for D-branes, one arrives at

(16)
for a D3-brane wrapped on the $`S^2`$. Using $`\tau _{D3}1/(\alpha ^2g_s)`$, one obtains

(17)
where $`\phi `$ has been rescaled to be a canonically normalized complex scalar field. We have derived (17) using a particular construction of light strings in four dimensions, but it is in fact typical for a variety of non-perturbative constructions, many of which are related to one another by string dualities.

Let us now consider the situation where after inflation, various moduli of string theory are rolling around, presumably before stabilizing in a preferred vacuum. The scalar $`\phi `$ in particular rolls around the complex plane, and if it comes close to the origin, then the light strings we have discussed become nearly tensionless. Based on our earlier discussion, we may expect appreciable string creation at such a moment. To estimate it, suppose $`\phi (t)\phi _0+\dot{\phi }t`$, where $`\phi _0`$ is the closest approach to the origin on the trajectory $`\phi (t)`$. Then

(18)
To estimate the amount of string creation, we want to find a complex solution $`t_{}`$ to $`\omega (t)^2=0`$. For highly excited strings, it is equivalent to solve $`\tau _{\mathrm{eff}}=0`$, which gives

(19)
Thus the convergence criterion on the total number of strings produced is

(20)
The minimum tension for the light strings is $`\tau _{\mathrm{min}}M|\phi _0|`$. We can estimate that $`\dot{\phi }M_{\mathrm{Pl},4}m_{\mathrm{inflaton}}`$: the $`\phi `$ oscillations are of amplitude $`M_{\mathrm{Pl},4}`$ and of frequency $`m_{\mathrm{inflaton}}`$. The condition (20) becomes

| $`1`$ | $`\begin{array}{c}<\hfill \\ \hfill \end{array}{\displaystyle \frac{\sqrt{M\phi _0^3}}{\dot{\phi }}}={\displaystyle \frac{\tau _{\mathrm{min}}^{3/2}}{M\dot{\phi }}}{\displaystyle \frac{\tau _{\mathrm{min}}^{3/2}}{MM_{\mathrm{Pl},4}m_{\mathrm{inflaton}}}}10^4\left({\displaystyle \frac{\tau _{\mathrm{min}}}{M_{\mathrm{Pl},4}^2}}\right)^{3/2},`$ |
| --- | --- |

(21)
or equivalently,

(22)
where we have assumed $`M10M_{\mathrm{Pl},4}`$ and $`m_{\mathrm{inflaton}}10^5M_{\mathrm{Pl},4}`$.

The estimates leading to (22) are inevitably fairly crude as long as we have not specified what the inflaton is and how much it overlaps with the scalar $`\phi `$ controlling the effective string tension. The result (22) is rather different from what we had for fundamental strings, and the reason is that we almost have a discontinuity in the first derivative of $`\omega (t)^2\tau _{\mathrm{eff}}M|\phi |M|\dot{\phi }t|`$.

Evidently, the non-perturbatively constructed strings do not have to dip very far below the Planck scale in order to violate this convergence criterion! So their participation in reheating is quite plausible. One would still have to answer the question of how dense in moduli space the loci of tensionless strings are. Quantifying this is technically difficult, but the fact that the loci are typically real codimension two is encouraging.

## 4 Backgrounds involving only Neveu-Schwarz fields

Let us now consider a background which might be more amenable to a worldsheet treatment.222This section represents work in progress with J. Friess and I. Mitra. As noted earlier, it is desirable to explore the possibilities for such vacua, since one would eventually want to go beyond the effective field theory method of estimating string creation that we described in section 2. Essentially all string theory compactifications give rise to the following terms in the effective action (written in four-dimensional string frame):

| $`S`$ | $`={\displaystyle \frac{1}{2\kappa _4^2}}{\displaystyle d^4x\sqrt{G}e^\varphi \left[R+(\varphi )^2\frac{1}{2}H_3^2\right]}`$ |
| --- | --- |
| $`e^\varphi `$ | $`={\displaystyle \frac{VolM_6}{(2\pi \sqrt{\alpha ^{}})^6}}{\displaystyle \frac{1}{g_s^2}}2\kappa _4^2=2\pi \alpha ^{}.`$ |

(23)
In four-dimensional Einstein frame, defined by $`ds_E^2=e^\varphi ds_{str}^2`$, we have instead

(24)
and the fundamental string tension as measured in Einstein frame is $`\tau _{\mathrm{eff}}=e^\varphi /2\pi \alpha ^{}`$. Except in type I string theory, $`H_3`$ is the four-dimensional part of the field strength of the Neveu-Schwarz two-form, and $`\chi `$ is defined via the equation $`d\chi =e^\varphi H_3`$.

The scalar $`\varphi `$ and the axion $`\chi `$ parametrize the upper half plane. As a first guess, we might suppose that there are FRW solutions where $`\chi +ie^\varphi `$ approximately describes a semi-circular geodesic in the upper half plane (with the usual metric), and that by making $`g_s`$ small, we can avoid significant back-reaction of the scalars on the geometry. Solving the scalar equations of motion in a background assumed to be flat Minkowski space, one finds

(25)
where $`\varphi _0`$ and $`t_0`$ are integration constants. The usual convergence criterion for the total number of strings created is $`t_0\sqrt{\tau _{\mathrm{eff}}(0)}\begin{array}{c}>\hfill \\ \hfill \end{array}1`$, or, more explicitly,

(26)
It would appear that by making $`\varphi _0`$ very negative (corresponding to extremely weak string coupling), with $`t_0`$ large compared to $`\sqrt{\alpha ^{}}`$ but fixed, we can get the bound in (26) to fail. This is a first indication that this type of solution is an interesting laboratory for string creation. It is nice that the wave equation for string modes is approximately a Mathieu equation, and that in the early and late time limits it is possible to construct infinite order adiabatic in and out vacua.

In fact, it is not difficult to promote the approximate solution described above to an exact solution of the equations of motion following from the two-derivative action:

| $`H_3`$ | $`=e^\varphi d\chi =hdx^1dx^2dx^3e^\varphi =e^{\varphi _0}\mathrm{sec}ht`$ |
| --- | --- |
| $`ds_{4E}^2`$ | $`=e^{\varphi _0}(\mathrm{sec}^2ht)\left({\displaystyle \frac{\mathrm{cos}\frac{ht}{2}+\mathrm{sin}\frac{ht}{2}}{\mathrm{cos}\frac{ht}{2}\mathrm{sin}\frac{ht}{2}}}\right)^\sqrt{3}dt^2+e^{\frac{\varphi _0}{3}}\left({\displaystyle \frac{\mathrm{cos}\frac{ht}{2}+\mathrm{sin}\frac{ht}{2}}{\mathrm{cos}\frac{ht}{2}\mathrm{sin}\frac{ht}{2}}}\right)^{1/\sqrt{3}}d\stackrel{}{x}^2.`$ |

(27)
Here $`t`$ runs from $`\pi /2h`$, where there is a curvature singularity, to $`\pi /2h`$, which is actually a non-singular infinite future when expressed in terms of proper time. This solution is not new: to my knowledge it first appeared in .

The main qualitative difference between the full solution (27) and the no-back-reaction approximation (25) is that the in vacuum is more difficult to define in (27), due to the singularity. Two possibly interesting approaches would be to replace the early time evolution by an appropriate instanton, or to find a CFT formulation in which $`\alpha ^{}`$ corrections are under control.

The solution (27) can be lifted to ten dimensions, and in the cases where it comes from IIA strings or $`E_8\times E_8`$ heterotic strings, there is a further lift to eleven dimensions, where the size of the eleventh dimension is controlled by the dilaton. When the lift is done to Horava-Witten theory, it describes two $`E_8`$ that start far apart at $`t=\pi /2h`$, move swiftly together to a minimum distance at $`t=0`$, and then move apart again. This is reminiscent of ekpyrosis, but actually it is quite different: in four-dimensional Einstein frame, there is never any contraction of the scale factor. There is also no singularity when the branes are at their distance of closest approach. The reason they don’t actually collide is that the four-form field $`G_4`$ is turned on in such a way as to make the $`E_8`$ branes repel. It is notable that in ekpyrosis as envisioned in , strings become tensionless at the instant that the branes collide. It seems extremely likely that string creation is important for understanding how the universe might pass through this singular event.

The point of involving only Neveu-Schwarz fields is that it is more likely that one can perform a worldsheet analysis of such backgrounds. In fact, the background we have considered resembles a limit of NS5-branes, doubly Wick rotated so that radius is replaced by time. The main stumbling block to finding a CFT description is that the geometry does not factorize in the way that the NS5-brane geometry does, into a three-sphere part, a linear dilaton background, and flat space.

## 5 Conclusions

Part of the appeal of the discussion in section 3 is that nearly tensionless strings arise in a variety of contexts, on loci in moduli space of real codimension two. Thus, independent of detailed ideas about how string theory connects to the Standard Model, it is plausible to have non-perturbatively constructed strings become temporarily light enough to play a role in post-inflationary cosmology. Because we remain ignorant of what string vacuum we might really be living in, it’s highly desirable to have some ideas about characteristic features of string theory which may become testable, which are difficult to contrive without invoking strings, and which are common to many backgrounds. Creation of light strings seems to qualify on all counts. On the other hand, the experimental constraints on precisely how the universe was reheated seem rather weak, so to really make the idea of post-inflationary string creation testable, one must gain a better idea of how such string creation would affect quantities that will be measured in the near future, like properties of the primordial spectrum of perturbations.

On the formal side, a better understanding of string creation from a worldsheet point of view is needed, and it would also be useful to learn more about the decay of highly excited strings. One reason that it is particularly important to have some examples of closed string creation described on the string worldsheet is that for bounds like (7) to be violated, one is usually forced to have string scale curvatures in the background.333The treatment of nearly tensionless strings in section 3 may be an exception to this. So a worldsheet treatment is needed to be sure one has $`\alpha ^{}`$ corrections under control, even before quantum creation processes are contemplated. I hope to report on some of these issues in the future.

## Acknowledgments

I thank the organizers of QTS3 for putting together a stimulating workshop. This work was supported in part by the Department of Energy under Grant No. DE-FG02-91ER40671, and also by a Sloan Fellowship.