text:
# Scaling in Ordered and Critical Random Boolean Networks

Permanent address: \]Physics Dept., Duke University, Durham, NC, 27514.

## Abstract

Random Boolean networks, originally invented as models of genetic regulatory networks, are simple models for a broad class of complex systems that show rich dynamical structures. From a biological perspective, the most interesting networks lie at or near a critical point in parameter space that divides “ordered” from “chaotic” attractor dynamics. In the “ordered” regime, we show rigorously that the average number of relevant nodes (the ones that determine the attractor dynamics) remains constant with increasing system size $`N`$. For critical networks, our analysis and numerical results show that the number of relevant nodes scales like $`N^{1/3}`$. Numerical experiments also show that the median number of attractors in critical networks grows faster than linearly with $`N`$. The calculations explain why the correct asymptotic scaling is observed only for very large $`N`$.

A random Boolean network (RBN) is a collection of $`N`$ binary logic gates, or nodes, wired together in a random fashion, with each node implementing a randomly chosen logical function of its inputs. RBNs are paradigms for systems in which excitatory and inhibitory interactions occur among a large set of interacting elements. One example of great current interest is the regulatory network that governs gene expression in a cell. It has been suggested that the distinct dynamical attractors of a single RBN be interpreted as distinct cell types carrying the same genetic information. Kauffman (1969) Surprisingly, RBN attractors can exhibit many features of biological cells, including stability against random external perturbations, qualitative change in response to special perturbations, and plausible scaling laws for numbers of attractors and attractor cycle lengths. Kauffman (1993) It therefore appears important to understand the behavior of RBNs as a first step in determining relevant global properties that might be probed in real gene expression experiments,

Even very simply constructed RBNs with deterministic updating rules can exhibit a rich set of dynamical behaviors. We focus here on the case in which each node has the same number of inputs, $`K`$. Fig. 1 shows an example with $`K=2`$. Each node $`i`$ implements a truth function $`F_i`$ (e.g. AND, XOR, etc. for K=2) that is chosen at random from a weighted distribution of all of the $`2^{2^K}`$ possible truth functions on $`K`$ binary inputs. On each discrete time step, the outputs are updated synchronously. Since the number of states of the system is finite (equal to $`2^N`$) and the system is deterministic, for any initial condition the network must eventually settle into a periodic attractor. We are interested in the behavior of large $`N`$ networks. How many attractors do they have? How many nodes typically participate in the attractor dynamics?

It is well known that tuning the probabilities of different $`F`$’s can produce an order-chaos transition. (See Aldana et al. (2002) for a thorough review.) In the ordered regime, almost all nodes are frozen and attractor cycles are short. In the chaotic regime, on the other hand, the number of fluctuating nodes is a finite fraction of $`N`$ and attractor cycles can be quite long. For large $`N`$, there is a narrow critical regime between these phases. The networks of greatest interest for biological systems are conjectured to lie near the critical regime, on the ordered side. Kauffman (1993)

In this Letter, we present numerical and analytic results that clarify the dynamical structure of RBNs in the ordered and critical regimes, revealing several surprising features: (1) for RBNs in the ordered regime the average number of relevant nodes (defined below) remains finite for $`N\mathrm{}`$ and they are organized into trivial loops; (2) in critical RBNs, the average number of fluctuating nodes grows like $`N^{2/3}`$; (3) the system sizes required to observe the asymptotic regime can be extremely large, especially for $`K=2`$; and (4) the median number of attractors in critical RBNs grows faster than linearly with $`N`$, at least for $`N`$ up to 1200. Both (2) and (4) contradict previous claims (Bastolla and Parisi (1998) and Bilke and Sjunnesson (2001), respectively) which we believe to have been based on studies that did not consider sufficiently large $`N`$. (4) also supersedes an old claim by one of us that the median number of attractors grows like $`\sqrt{N}`$ in critical networks. Kauffman (1969)

The concepts of “relevant” nodes Flyvbjerg and Kjaer (1988) and “canalizing inputs” Kauffman (1969) are essential to our analysis. In any given network, there may be nodes whose outputs are frozen at the same value on every attractor. Such nodes serve only to fix inputs to other nodes and are otherwise “irrelevant”. There may also be nodes whose outputs go only to irrelevant nodes. These are also classified as irrelevant. Though they may fluctuate, they act merely as slaves to the nodes that determine the attractor cycle.

Almost all the irrelevant nodes can be found as follows. Bilke and Sjunnesson (2001) One first identifies “fixed” nodes whose outputs are entirely independent of their inputs. One then uses an iterative procedure to identify nodes that must be frozen because their outputs depend only on inputs from other frozen nodes. We call all frozen nodes identified this way “clamped”, let $`s`$ denote their number, and define $`uNs`$. A similar iterative procedure is then used to remove (or “prune”) nodes with no relevant outputs. For example, in Fig. 1, even if nodes $`15`$ are all unclamped, nodes $`68`$ can be pruned. For the purposes of this paper we designate all nodes that are neither clamped nor pruned via the described procedure as “relevant”, and let $`r`$ denote their number. (Additional nodes may be frozen due to correlations between two or more unclamped inputs Bilke and Sjunnesson (2001), so $`r`$ is greater than or equal to the number of truly relevant nodes.)

A canalizing input to a given node is one that can be set to a value that determines the output, independent of the other inputs. (For example, if either input to an OR function is set to 1, the output is determined.) Following Flyvbjerg and Kjaer (1988), we define a set of parameters $`p_k`$ as follows. For a randomly selected $`F`$, fix a randomly selected $`Kk`$ inputs at arbitrarily chosen values. $`p_k`$ is the probability that those input values are collectively canalizing; i.e., that $`F`$ is independent of the $`k`$ remaining inputs. Note that $`p_0=1`$ (fixing all the inputs certainly determines the output) and $`p_K`$ is the probability that a node is fixed (i.e., that its output is independent of all of its inputs).

The order-chaos transition can be observed by tuning the $`p_k`$’s. One simple way to do this is to assign to each $`F`$ a probability that depends only on the number of $`1`$ in the output column of its truth table. For $`p[0,1]`$, we let the probability that a node has truth function $`F`$ be $`p^{k_1}q^{Kk_1}`$, where $`q=(1p)`$. It is straightforward to check that this parametrization corresponds to $`p_k=p^{2^k}+q^{2^k}`$. With this weighting of the $`F`$’s, the transition occurs at $`2Kpq=1`$, with $`2Kpq<1`$ corresponding to the ordered regime Derrida and Flyvbjerg (1987); Derrida and Pomeau (1986); Derrida and Weisbuch (1986). For our numerical studies, we set $`K=2`$ and vary $`p`$.

We have carried out two types of numerical experiment on $`K=2`$ networks. First, for 1000 networks at each of five $`p`$ values, we determine $`u`$ and $`r`$. Fig. 2(a) shows the measured $`r`$ for $`N`$ up to $`3\times 10^6`$. It appears that $`r`$ approaches a constant at large $`N`$ for $`p`$ in the ordered regime. At the critical value $`p=1/2`$, the data are inconclusive, showing significant curvature on the log-log plot out to the largest $`N`$ we have studied. They are consistent, however, with an asymptotic scaling law $`rN^{1/3}`$. Panel (b) shows the $`u`$ for $`p=1/2`$, indicating a clearer power-law scaling $`uN^{2/3}`$.

Second, for at least 1000 networks at each $`p`$, we attempt to measure the number of attractors, $`A`$, on the set of relevant nodes. In some networks, however, $`A`$ is prohibitively large, making measurements of $`A`$ difficult. It is much easier to measure the median, $`\widehat{A}`$, since one need not continue to count attractors in a given network after the count has exceeded the median. To count attractors, we repeatedly choose random initial conditions and identify the attractor reached. If 1000 consecutive attempts yield no new attractor, we record the number of attractors found and move on to another network. This gives a lower bound on $`A`$ for each network, and hence a lower bound on $`\widehat{A}`$.

Fig. 3 shows the results for $`N`$ up to 1200. We note that in Bilke and Sjunnesson (2001), where measurements of $`A`$ were sought, it was not possible to consider nets larger than $`N=144`$. From Fig. 3, however, it is clear that any extrapolation based on data for $`N`$ smaller than about 500 at the critical point is suspect. We see a faster than linear rise in the median $`A`$ above $`N500`$, which almost certainly implies a faster than linear rise in $`A`$ as well. Moreover, Fig. 2(a) strongly suggests that one must study $`N>10^6`$ to observe the true asymptotic behavior!

A simple calculation provides a rigorous upper bound on $`u`$ and explains why asymptotic scaling sets in only for very large $`N`$. A technical note: in the networks studied numerically, gates were not permitted to have two inputs from the same node. The analysis below ignores this constraint, which yields $`𝒪(1/N)`$ effects. In all cases, self inputs are allowed.

Let $`P(u)`$ be the probability that a randomly selected network has $`u`$ unclamped nodes. To compute $`P(u)`$, we should count the networks with $`u`$ unclamped nodes and divide by the total number of networks $`T=[\tau N^K]^N`$, where $`\tau =2^{2^K}`$ is the number of possible truth functions for each node. We begin by considering a quantity $`P^{}(u)`$ that is guaranteed to exceed $`P(u)`$:

$`P^{}(u)={\displaystyle \frac{1}{T}}`$ $`C(N,u)`$ $`\left[{\displaystyle \underset{k=0}{\overset{K}{}}}C(K,k)u^ks^{Kk}\tau p_k\right]^s`$
$`\times `$ $`\left[{\displaystyle \underset{k=0}{\overset{K}{}}}C(K,k)u^ks^{Kk}\tau q_k\right]^u,`$ (1)
where $`q_k1p_k`$ and $`C(m,n)`$ is the number of combinations of $`n`$ objects drawn from $`m`$. The factor $`C(N,u)`$ is the number of ways of having $`u`$ unclamped nodes. The first sum counts the ways that a node can be clamped, weighted by the probability of its truth function: the node in question can have up to $`k`$ unstable inputs as long as the other $`Kk`$ are collectively canalizing, which occurs with probability $`p_k`$. Simlarly, the second sum counts the ways a node can be unclamped.

$`P^{}(u)`$ overcounts the probability of having $`u`$ unclamped nodes because it ignores the constraint that all clamped nodes must be traceable through a sequence of inputs back to a fixed node. That is, the first sum overcounts the number of ways that $`s`$ clamped nodes can be wired, as it includes graphs in which a subset of the $`s`$ nodes collectively clamp each other without any connection to a fixed node. For example, suppose that in Fig. 1 nodes $`1`$ and $`4`$ are fixed, node $`9`$ implements OR, and the output of $`10`$ is simply equal to its input from $`9`$. Nodes $`2`$, $`3`$, and $`58`$, are clearly clamped through inputs that can be traced back to $`1`$ and $`4`$. The sum then counts this network as a possible arrangement of the clamped nodes for $`u=0`$ because it is self-consistent to assume that both $`9`$ and $`10`$ are clamped. (Note that after a few time steps 9 and 10 will either both be stuck on 0 or both on 1.) However the iterative procedure defined above would (correctly) not identify $`9`$ and $`10`$ as clamped, so its inclusion in $`P^{}(0)`$ constitutes overcounting. Note that the same network is also (properly) counted in $`P^{}(2)`$. Thus $`P^{}(u)>P(u)`$ for all $`u`$, implying that $`_ug(u)P^{}(u)_ug(u)P(u)`$ for any $`g(u)0`$.

If $`u`$ is small compared to $`N`$, a useful approximation to Eq. (Scaling in Ordered and Critical Random Boolean Networks) can be obtained. Using $`s=Nu`$ and using Stirling’s formula to simplify the binomial coefficients, together with the identity $`(1+x/n)^n=\mathrm{exp}[x(1(x/n)/2+(x/n)^2/3+\mathrm{})]`$, we find

$$P^{}(u)\frac{1}{\sqrt{2\pi u(1u/N)}}\mathrm{exp}\left[\theta _1u\theta _2\frac{u^2}{N}\theta _3\frac{u^3}{N^2}\right],$$
(2)
where

$`\theta _1`$ $`=`$ $`Q1\mathrm{ln}Q,`$ (3)
$`\theta _2`$ $`=`$ $`(1Q)(1+Q2K+2Q_2)/2,`$ (4)
$`\theta _3`$ $`=`$ $`\theta _2P_3Q_3+(KQ)P_2+{\displaystyle \frac{1}{2}}Q_2^2`$
$`[K1+2(KQ)^3]/6,`$ (5)
with $`QKq_1`$, $`Q_nC(K,n)q_n/Q`$, and $`P_nC(K,n)p_n`$. The calculation involves only straightforward algebra and the assumption that terms of order $`u^4/N^3`$ can be neglected in the exponent, which is justified whenever $`\theta _1`$, $`\theta _2`$, and $`\theta _3`$ are all positive and at least one of them is nonzero. Under such circumstances, and in the limit of large $`N`$, $`P^{}(u)`$ is exponentially strongly suppressed for $`u`$’s larger than order $`N^{2/3}`$, whereas the neglected terms would only become relevant for $`u`$’s of order $`N^{3/4}`$.

$`Q`$ is equivalent to the order parameter defined by Flyvbjerg in Flyvbjerg and Kjaer (1988), where it was also argued that $`Q=1`$ marks the critical boundary. Eq. (2) proves that for all $`Q<1`$ the $`P(u)`$ decays exponentially with a decay length independent of $`N`$. It also strongly suggests that this is not the case at $`Q=1`$, though the computation only gives an upper bound. The fact that $`\theta _1`$ and $`\theta _2`$ both vanish at $`Q=1`$ is a surprising result that affects the scaling of $`r`$ in critical networks, as we shall see below.

In the ordered regime $`Q<1`$, we have $`\theta _1>0`$, so higher order terms are irrelevant at sufficiently large $`N`$ and $`P^{}(u)`$ becomes independent of $`N`$. Thus $`u`$ for asymptotically large $`N`$ is bounded above by a constant. For the case where $`p_k`$ is determined simply by the one parameter $`p`$, we have $`q_1=2p(1p)`$ and the critical $`p`$ satisfies $`2Kp_c(1p_c)=1`$, consistent with previous studies of the effect of $`p`$. For $`p`$ near $`p_c`$ we get

$$\theta _1\{\begin{array}{cc}8(pp_c)^4\hfill & \mathrm{for}K=2\hfill \\ 2K(K2)(pp_c)^2\hfill & \mathrm{for}K>2.\hfill \end{array}$$

For $`K=2`$, then, $`\theta _1`$ is quite small even for $`p`$ relatively far from $`p_c=1/2`$. For example, $`p=0.55`$ gives $`\theta _15\times 10^5`$. The asymptotic behavior of the system should only be apparent when $`u`$ is bigger than $`1/\theta _1`$ and hence when $`N`$ is substantially larger than that. The vertical dashed lines in Fig. 2(a) are drawn at $`N=3/\theta _1`$, which is seen to give an accurate indication of where effects associated with the ordered regime become important.

Moreover, for very large $`N`$, we may obtain an accurate approximation to $`P(u)`$ simply by normalizing $`P^{}(u)`$, since the overcounting factor relating the two approaches some nonzero constant at $`u=0`$. This allows us to compute $`r`$ (not just an upper bound on it) as follows. Consider the unclamped nodes and the links between them. These form a random graph subject only to the constraint that each node has between 1 and $`K`$ inputs. But from the second sum in Eq. (Scaling in Ordered and Critical Random Boolean Networks) we know the relative probabilities that an unclamped node will have $`k`$ unclamped inputs. Specializing to $`K=2`$ for simplicity, the average number of inputs per node for fixed $`u`$ is

$`z`$ $`=`$ $`\left(Qsu+2q_2u^2\right)/\left(Qsu+q_2u^2\right)`$
$``$ $`1+{\displaystyle \frac{q_2}{Q}}{\displaystyle \frac{u}{N}}+\left({\displaystyle \frac{q_2}{Q}}{\displaystyle \frac{q_2^2}{Q^2}}\right){\displaystyle \frac{u^2}{N^2}}.`$ (6)
Now in the ordered regime, where the exponential cutoff in $`P(u)`$ is independent of $`N`$, we have $`𝑑uP(u)u/N0`$ for all $`u`$ as $`N\mathrm{}`$ and hence $`z1`$. Thus the network of unclamped nodes at large $`N`$ is just a randomly wired $`K=1`$ network with no fixed nodes. Exact combinatorics for $`K=1`$ random graphs show that the expected number of loops of size $`n`$ in a network of size $`u`$ is $`L_u(n)=u!/[n(un)!(u1)^n].`$ Flyvbjerg (1988) Since all nodes not in loops can be pruned, we obtain

$$r=\underset{u=0}{\overset{N}{}}\underset{n=1}{\overset{u}{}}nL_u(n)P(u)0.7/\sqrt{\theta _1}.$$
(7)
The horizontal dashed lines on the far right in Fig. 2(a) mark this prediction, which agrees well with the data.

For the critical case $`Q=1`$, we have

$$P^{}(u)[2\pi u(1u/N)]^{1/2}\mathrm{exp}\left[\theta _3u^3/N^2\right],$$
(8)
which implies that $`u`$ cannot grow faster than $`N^{2/3}`$. This contradicts a previous argument suggesting $`uN^{3/4}`$ Bastolla and Parisi (1998). Note however, that the corrections from terms of order $`u^4/N^3`$ cannot be neglected unless $`(N^{2/3})^4/N^3=N^{1/9}<<1`$, suggesting that the critical scaling emerges only for $`N>10^9`$, as confirmed by Fig. 2.

Naive normalization of $`P^{}(u)`$ to get $`P(u)`$ yields $`uN^{2/3}`$ and Fig. 2(b) indicates that this is correct. Given the scaling law and the fact that paths connecting relevant nodes cannot lead to dead ends, the theory of directed random graphs Newman (2001) and Eq. (6) can be used to show that $`rN^{1/3}`$, consistent with the numerical trend seen in Fig. 2(a). Unfortunately, naive normalization does not yield an accurate prediction for the full form of $`P(u)`$ because the $`u`$ dependence of the factor $`f(u)P(u)/P^{}(u)`$ becomes important. Determination of the precise form of $`f(u)`$ is beyond the scope of this work.

We now address the question of the number of attractors. For the case $`Q<1`$ and very large $`N`$, the relevant nodes form trivial loops with only two possible truth functions: $`F(\sigma )=\sigma `$ (“identity”) and $`F(\sigma )=1\sigma `$ (“not”). For $`n`$ prime, a loop of size $`n`$ has either $`(2^n2)/n`$ or $`(2^n2)/2n`$ attractors, depending on whether the number of nots in the loop is even or odd. For $`n`$ not prime, the number of attractors is slightly larger. Flyvbjerg (1988) Now $`L_u(n)(2^n2)/n`$ is extremely sharply peaked very close to $`n=u/2`$. $`A_u`$ is thus dominated by rare networks that contain a relevant loop of size $`u/2`$, which gives $`A_u2^{cu/2}`$, with $`c`$ a constant of order unity. Naive averaging over $`P(u)`$ gives a divergent answer; a consistent calculation would require inclusion of the terms in $`P(u)`$ that cause rapid decay for $`u`$ of order $`N`$.

The median number of attractors, $`\widehat{A}`$, can be estimated as $`2^{\widehat{b}}/\widehat{b}`$, where $`b`$ is the size of the biggest loop of relevant nodes in a given network. For fixed $`u`$, the probability of occurence of $`b`$ is $`P_u(b)=L_u(b)_{n=b+1}^u\left[1L_u(n)\right]`$. To find $`\widehat{b}`$, we numerically solve $`_{b=0}^{\widehat{b}}_{u=0}^{\mathrm{}}P(u)P_u(b)=1/2`$. The results for the networks studied in Fig. 3 are $`\widehat{A}=3`$, $`5`$, and $`14`$ for $`p=0.7`$, $`0.65`$, and $`0.6`$, respectively. This result and its rough agreement with the data (see Fig. 3) confirms that $`\widehat{A}`$ approaches a constant at large $`N`$ for any sub-critical $`p`$ lyn  and illustrates the qualitative distinction between $`A`$ and $`\widehat{A}`$.

In closing, we note that the classification of nodes as relevant is robust across a large class of updating rules. The results on $`r`$ reported here apply to asynchronous models and to models in which there is a stochastic time delay in the updating of each node, as long as the update, whenever it occurs, is always accurate. Finally, we note that the numbers of genes in real cells are on the order of $`10^4`$, which our results show may be too small to exhibit asymptotic large $`N`$ behavior. Thus, networks with canalization parameters that nominally place them in the ordered regime will exhibit features of critical networks, which may be important for biological function.

Acknowledgments: We thank M. Newman, M. Girvan, S. Coppersmith, R. Palmer, and H. Flyvbjerg for helpful communications. J.E.S.S. thanks Bios Group, the Santa Fe Institute, and the Smith Faculty Enrichment Fund at Duke University for financial and facilities support.


format: text

text:
# Qubit versus bit for measuring an integral of a classical field

## Abstract

Methods for measuring an integral of a classical field via local interaction of classical bits or local interaction of qubits passing through the field one at a time are analyzed. A quantum method, which has an exponentially better precision than any classical method we could see, is described.

Although the retrievable information content of a number of qubits is essentially equal to the information content which can be stored in the same number of bits Ho , qubits are more efficient for many specifically tailored tasks. Recently, Galv$`\stackrel{~}{\mathrm{a}}`$o and Hardy GH  (GH) found one such task: pinpointing a particular property of an integral of a classical field is possible using a single qubit, but impossible with a classical bit, or even with many bits. This motivated the current work in which we found that qubits are more efficient than bits for a more general task, namely, the measurement of the integral itself.

The task is to measure the integral of a classical field,

$$I=_A^B\varphi (x)𝑑x,$$
(1)
via local interactions along the path from $`A`$ to $`B`$. In one case $`N`$ bits go one at a time through the path interacting with field, and in the other case, $`N`$ qubits pass through the field instead. In this task there are no constrains on the complexity of local interactions between the bits (qubits) and the field. The limitations which we investigate are due to the complexity of the carrier of the information about the field which goes along the path.

In order to make the comparison between quantum and classical methods simpler, we consider a non-negative classical field and assume that it is known that the order of magnitude of $`I`$ is $`M`$. The classical method is as follows. The bit starts at $`A`$ in the state 0 and it flips with a probability proportional to the strength of the field. Once flipped to 1, it cannot flip back. The probability of a flip along the path is

$$p=1e^{\lambda I},$$
(2)
where $`\lambda `$ has to be optimized for getting the best precision.foot1  The exact optimization depends on our a priory information about the probability for different values of $`I`$ and the choice of the particular aspect of the precision of the measurement we wish to optimize. The precision of measuring $`p`$ achieved after sending $`N`$ bits can be estimated as $`\mathrm{\Delta }p=\frac{\sqrt{p(1p)}}{\sqrt{N}}`$. Thus, essentially for all reasonable approaches, the uncertainty of measuring $`I`$ is of the order of $`\frac{M}{\sqrt{N}}`$.foot2

In the quantum method, we arrange an interaction which leads to a qubit “rotation” proportional to the strength of the field, for example, a spin precession in a magnetic field. The method works equally well for negative and general classical fields. The only important parameter is the range of possible values of the integral. The strength of the coupling should be such that for every value of the integral we will get different final state of the qubit. Thus, we have to choose the strength of the coupling to the field such that a rotation of more than $`2\pi `$ will not be probable.

When we send $`N`$ spin-$`\frac{1}{2}`$ particles through the field, they all rotate by the same angle. If we start with all spins pointing in, say $`x`$ direction, and the magnetic field in the $`z`$ axis, the direction of the spin in the $`xy`$ plane will yield the value of the integral.

The precision of measuring direction of $`N`$ parallel spins is proportional to $`N^{\frac{1}{2}}`$. This is the same dependence as in the probabilistic methods with bits. However, Peres and Scudo PS  (see also Bagan et al. Ba ) showed that taking $`N`$ entangled spins one can reduce the uncertainty to be proportional to $`N^1`$. They optimized defining the direction in three dimensions with $`N`$ spin-$`\frac{1}{2}`$ particles, while in the present problem we need to find direction only in two dimensions. So, some further optimization is possible, but the uncertainty remains proportional to $`N^1`$. This concludes the description of our quantum method I for measurement of integral $`I`$.

A crucial requirement for the advantage of qubits versus bits is that we are only allowed to send bits and qubits one at a time. If we send all bits together, we can build a probabilistic counter which can register up to $`2^N`$ counts. The counter makes counts while moving in the field with the probability for a count proportional to the field. The uncertainty in the total number of counts is of the order of the square root of this number. Therefore, adjusting the strength of the interaction such that the expectation value of the number of counts is of the order of $`2^N`$, we can get the uncertainty in the measurement of $`I`$ proportional to $`2^{\frac{N}{2}}`$. The uncertainty is exponentially smaller than in the quantum method described above.

Of course, if qubits are allowed to pass along the field together, the precision that can be achieved in the measurement of the integral is not smaller than that in the classical case. But can exponentially small uncertainty be achieved with qubits passing through the field one at a time? As we will show below, the answer is yes, if we allow various strengths of the interaction between the qubits and the field.

On the other hand, we do not see how this freedom might lead to a significant advantage in the classical case. Indeed, it is intuitively clear (especially after the GH proof for a particular case) that a bit passing through the field cannot have a deterministic information about the integral of the field. Then, the only way we can imagine for storing information in the bit is in the value of its probability to be 1. When we are given $`N`$ bits, the uncertainty in the measurement of the probability decreases as $`N^{\frac{1}{2}}`$. It seems to us highly implausible that there is a classical method which do better than this.

Now we turn to the description of the method (referred hereafter as quantum method II) which employs qubits interacting with various strengths with the field. The strengths depend only on the number of the qubit passing through the field and are fixed before the experiment. No additional structure is required for the qubit.

The basis of this method is the result of Galv$`\stackrel{~}{\mathrm{a}}`$o and Hardy GH . They considered a particular case in which

$$I=m\alpha ,$$
(3)
where $`\alpha `$ is known and $`m`$ is an integer. GH found a method which allows with a single qubit to answer the question: is $`m`$ even or odd?foot3  They achieved the goal by tuning the strength of the interaction in such a way that $`I=\alpha `$ yields a rotation by $`\pi `$. Thus, for an odd integer the spin flips, and for an even integer, it returns to its initial state.

If we send a number of qubits, one after the other, we can modify this procedure to find $`m`$ itself. To this end it is arranged that the qubits we send interact with the field with different strengths: the first as in the GH protocol, the second with half of the strength of the first, the next with the half of the strength of the preceding, etc. In the first step we find the last digit of $`m`$ written in the binary way. In the next step we find the preceding digit, and thus, after $`k`$’th step we find $`m(\mathrm{mod2}^k)`$.

The procedure works in the following way. If the last digit is zero, then in the second step we repeat the GH protocol with half of the strength of the interaction. Since now we know that $`I=m^{}2\alpha `$, the protocol determines the last digit of $`m^{}`$ which is second digit from the end of $`m`$. If the last digit of $`m`$ is $`1`$, in the second step we should modify the procedure by additional rotation of the spin by the angle $`\theta _2=\frac{\pi }{2}`$. In the $`k`$’th step we should compensate for all non-zero digits, as follows:

$$\theta _k=\underset{i=1}{\overset{k1}{}}\frac{\pi d(i)}{2^{(ki)}}.$$
(4)
where $`d(i)`$ is the value of digit number $`i`$ from the end. The method yields one digit for each qubit and yields zeros once the whole number is written.

In a general case, when our only prior information about $`I`$ is its order of magnitude, we can combine the two methods we described above. We chose $`\alpha `$ and, using method I, we find the reminder $`\beta `$ of the division of $`I`$ by $`\alpha `$. Then, using method II, we find $`m`$:

$$I=m\alpha +\beta .$$
(5)
After measuring $`\beta `$ with a good precision we apply method II for finding $`m`$. The reminder $`\beta `$ requires additional correction angles in the application of the second method. For the first step, the correction angle is $`\theta _1=\frac{\pi \beta }{\alpha }`$ and, in general, for $`k`$’th step, the correction angle is:

$$\theta _k=\underset{i=1}{\overset{k1}{}}\frac{\pi d(i)}{2^{(ki)}}\frac{\pi \beta }{2^{k1}\alpha }.$$
(6)
The requirement for choosing $`\alpha `$ is that we will have enough qubits to find all digits of $`m`$. We get high probability for that if

$$\alpha =\frac{10M}{2^{NN_0}},$$
(7)
where $`N_0`$ is the number of qubits used in the measurement of $`\beta `$. If the probability of error in the measurement of $`m`$ is negligible (when $`\beta `$ is measured with the high precision), then the uncertainty in the value of $`I`$ is, essentially, the uncertainty in the measurement of $`\beta `$ which, in the best case (Peres-Scudo method), is of the order of $`\frac{\alpha }{N_0}`$.

However, as the anonymous referee pointed out, the error in the measurement of $`m`$ turns out to be not too large even if there is a large error in the measurement of $`\beta `$. In fact, it is more effective not to “waste” qubits on measurement of $`\beta `$. The quantum method II works well by itself even in a general case. Now $`\alpha `$ becomes much smaller:

$$\alpha =\frac{10M}{2^N}.$$
(8)
In order to estimate the precision of the quantum method we calculate probability of a particular reading $`m`$ given the actual value of $`I`$. We will show that it is rapidly decreasing function of the difference between the actual value of the integral $`I`$ and the read out of the device $`\stackrel{~}{I}=m\alpha `$. In our procedure, the last digit of $`m`$ is specified by the spin measurement. The spin rotates by the angle $`\mathrm{\Theta }_1=\frac{I\pi }{\alpha }`$ and is found in the direction specified by the angle $`\stackrel{~}{\mathrm{\Theta }}_1=\frac{\stackrel{~}{I}\pi }{\alpha }`$. Thus, the probability for this particular outcome of the spin measurement which specifies the last digit is $`p_1=\mathrm{cos}^2\frac{\mathrm{\Theta }_1\stackrel{~}{\mathrm{\Theta }}_1}{2}=\mathrm{cos}^2\frac{(I\stackrel{~}{I})\pi }{2\alpha }`$. The digit $`k`$ from the end, $`k>1`$, is specified by the measurement of the spin which is rotated by the angle

$$\mathrm{\Theta }_k=\frac{I\pi }{2^{k1}\alpha }+\theta _k$$
(9)
and found in the direction specified by the angle

$$\stackrel{~}{\mathrm{\Theta }}_k=\frac{\stackrel{~}{I}\pi }{2^{k1}\alpha }+\theta _k.$$
(10)
Therefore, the probability for this particular outcome of the spin measurement which specifies the $`k`$’th digit is $`p_k=\mathrm{cos}^2\frac{(I\stackrel{~}{I})\pi }{2^k\alpha }`$. Thus, the probability for readout $`\stackrel{~}{I}`$ given the actual value of the integral $`I`$ is:

$$p(\stackrel{~}{I}|I)=\underset{k=1}{\overset{N}{}}\mathrm{cos}^2\frac{(I\stackrel{~}{I})\pi }{2^k\alpha }.$$
(11)
Denoting the error $`\delta I`$ and substituting the value of $`\alpha `$ from (8), the probability of the error is

$$p(\delta I)=\underset{k=1}{\overset{N}{}}\mathrm{cos}^2\frac{\delta I\pi }{2^{kN}10M}.$$
(12)
This function vansihes for $`\delta I=n\alpha `$ for integer $`n`$, except for $`n=0`$, where it has maxima which is equal 1. It has to be the so, because the method yields no errors for $`I=m\alpha `$. The function has local maxima at $`\delta I=(n+\frac{1}{2})\alpha `$, except for $`n=1,0`$.

The read-out value of the integral, $`\stackrel{~}{I}`$, might have only discrete values, $`m\alpha `$. Therefore, the error of the order of $`\alpha `$ is unavoidable. In the worse case, $`I\mathrm{mod}\alpha =\frac{\alpha }{2}`$, but even in this case we have only small probability to get the error which is an order of magnitude larger than $`\alpha `$. For large $`N`$, this probability has almost no depndence on $`N`$; we calculate it for $`N=30`$:

$`p(\delta I>10\alpha )=1p(\delta I<10\alpha )`$
$`=1{\displaystyle \underset{n=9}{\overset{10}{}}}{\displaystyle \underset{k=1}{\overset{30}{}}}\mathrm{cos}^2{\displaystyle \frac{(\frac{1}{2}+n)\alpha \pi }{2^{kN}10M}}`$
$`=1{\displaystyle \underset{n=9}{\overset{10}{}}}{\displaystyle \underset{k=1}{\overset{30}{}}}\mathrm{cos}^2{\displaystyle \frac{(\frac{1}{2}+n)\pi }{2^k}}0.019.`$ (13)
Thus, we should expect an error of the order of $`\alpha `$ (it is of the order of $`10^7`$ in our case) and not significantly larger.

In order to illustrate our result we performed computer simulation of classical ($`30`$ bits) and quantum ($`30`$ qubits) measurements of the integral for ten different fields. We took $`I_n=(n\pi )\mathrm{mod10}`$, $`n=1,2,\mathrm{}10`$. We assumed that the order of magnitude of the integral is given, $`M=5`$, and we choose parameter $`\lambda `$ of the classical method such that the precision of the measurement of $`I`$ for $`I=M`$ is optimized. The uncertainty in classical measurement can be estimated as

$$\mathrm{\Delta }I=\frac{\sqrt{e^{\lambda I}1}}{\lambda \sqrt{N}}.$$
(14)
For $`I=M`$ it has a minimum around $`\lambda =\frac{1.2}{M}`$ and for parameter we chose, the uncertainty $`\mathrm{\Delta }I`$ is of the order of 1.

The results are shown in Table 1. We see that the error in classical method is indeed of the order of 1 and the quantum error is of the order of $`10^7`$.

The technology today is far from getting this exponential advantage of the quantum method which requires stability and high precision at very large range of the interaction strength. Also, preparation of initial entangled states and complicated collective measurements of Peres-Scudo measurements are difficult for experimental implementation. However, recent progress in quantum information experiments allows us to believe in the prospects of at least partial implementation of our proposal which will manifest advantage of the quantum method.

It is a pleasure to thank Yakir Aharonov, Shmuel Nussinov, Benni Reznik, Boris Tsirelson, and Stephen Wiesner for helpful discussions. This research was supported in part by grant 62/01 of the Israel Science Foundation and by the Israel MOD Research and Technology Unit.


text:
# A Method for Measuring Emittance in 𝑒⁺𝑒⁻ Colliding Beams

(4 December 2002)

## Abstract

We have developed techniques that allow simultaneous measurement of the spatial size of the luminous colliding beam region and the angular spread of beams in collision using $`e^+e^{}\mu ^+\mu ^{}`$ events. These are demonstrated at the CLEO interaction point of the Cornell Electron-Positron Storage Ring, CESR, taking advantage of the small and well understood resolution of the CLEO tracking system. These measurements are then used to extract the horizontal beta, horizontal emittance and the vertical emittance and search for dynamic effects at CESR.

preprint: CLNS 02/1812preprint: WSU-HEP-0210

The measurement of the emittance of colliding beams is a difficult and important problem as many methods are destructive or disruptive to the beams or involve extrapolation using a theoretical accelerator transport matrix seeman . Here we present a method that relies on the tracking system of a high energy physics experiment done while the beams are producing luminosity without disruption or use of a transport matrix. The basic technique is to select an ensemble of $`e^+e^{}\mu ^+\mu ^{}`$ events with which we simultaneously measure the size of the luminous region and the intrinsic angular spread of the outgoing particles. At CESR and CLEO the resolution is precise enough that we can extract the the horizontal beta at the interaction point, $`\beta _x^{}`$, the horizontal emittance, $`ϵ_x`$, and, combined with our previously described technique for measuring vertical beta at the interaction point, $`\beta _y^{}`$ hourglass , the vertical emittance, $`ϵ_y`$.

Assuming that dispersion and coupling are negligible, the beam parameters $`\beta _x^{}`$, $`\beta _y^{}`$, $`ϵ_x`$, and $`ϵ_y`$ are related to the physical observables $`\sigma `$, the Gaussian width of the beam at the interaction point, and $`\sigma ^{}`$, the Gaussian width of the angular spread of the beam at the interaction point, by the following equations:

$`\sigma _x`$ $`=`$ $`\sqrt{\beta _x^{}ϵ_x}`$ (1)
$`\sigma _y`$ $`=`$ $`\sqrt{\beta _y^{}ϵ_y}`$ (2)
$`\sigma _x^{}`$ $`=`$ $`\sqrt{ϵ_x/\beta _x^{}}`$ (3)
$`\sigma _y^{}`$ $`=`$ $`\sqrt{ϵ_y/\beta _y^{}}.`$ (4)
A simultaneous measurement of both $`\sigma `$ and $`\sigma ^{}`$ in the collision region would yield measurements of $`\beta ^{}`$ and $`ϵ`$ via the obvious algebraic manipulation.

Our previous work has shown how the size of the luminous region can be measured using $`e^+e^{}\mu ^+\mu ^{}`$ events hourglass . Briefly a fiducial box is centered on the measured center of the luminous region. The average positions of tracks passing through the box are used to measure the size and shape of the luminous region. Tracks that pass through the box nearly perpendicular to one of the sides are most useful for making a measurement. At the CESR-CLEO interaction point the luminous region is roughly 7 $`\mu `$m high and 300 $`\mu `$m wide. Thus we select more useful tracks for the horizontal measurement than the vertical. We have previously used this technique to measure $`\beta _y^{}`$ and observe the hourglass effect at the CESR-CLEO interaction point hourglass .

We measure the angle between the projected $`\mu ^+`$ and $`\mu ^{}`$ tracks by considering the sum of the momenta of the two tracks divided by the beam energy to obtain a vector. If the two tracks are coming from the same point and have nearly the same momentum as the beam this vector is equivalent to the minor arc between the colliding tracks, given by $`E_{\mathrm{beam}}\mathrm{sin}\theta `$, where $`\theta `$ is the central angle between the colliding tracks. This angle is small and we approximate $`\mathrm{sin}\theta `$ by $`\theta `$. We then define the angle between the two tracks in the horizontal projection as the horizontal momentum sum divided by beam energy and similarly in the vertical direction. Schematically this is shown in Figure 1.

A crossing angle between the beams causes the distribution of the angle to have a mean different from zero, an angular spread of the beam, $`\sigma ^{}`$, gives an intrinsic underlying width, and detector momentum resolution broadens the distribution.

CESR has been described in detail elsewhere CESR . All the data used in this measurement were taken at an $`e^+e^{}`$ collision energy of 10.58 GeV with bunch currents in the range of 2.5 to 7.0 mA over a four month period in late 1998 and early 1999. The CLEO detector has also been described in detail elsewhere CLEO . All of the data used in this measurement are taken in the CLEO II.V configuration which includes a silicon strip vertex detector (SVX) which is crucial to the measurement of the size of the luminous region. This detector consists of three layers of silicon wafers arrayed in an octagonal geometry around the interaction point. The first measurement layer is at a radius of 2.3 cm and the wafers are read out on both sides by strips which are perpendicular to each other. The readout strips have a pitch of about 100 $`\mu `$m; with charge sharing, the detector has an intrinsic per-point resolution of better than 20$`\mu `$m both longitudinally and in the plane transverse to the beam direction. The rest of the tracking system consists of a small 10 layer drift chamber, and a large 51 layer drift chamber which has an outer radius of 1.5 m. The tracking system is in a 1.5 Tesla magnetic field provided by a superconducting solenoid.

$`e^+e^{}\mu ^+\mu ^{}`$ events are easily selected in CLEO by choosing events with two and only two tracks each with momentum near the beam energy and a small energy deposit in the electromagnetic calorimeter. We chose tracks with 20 or more hits in the main drift chamber, and at least two silicon vertex detector hits in the transverse and longitudinal views. We require that the tracks have opposite charge and that those used for the measurement of the luminous region have at least three silicon vertex detector hits in one of the two views. Each point for measuring the $`\sigma ^{}`$ is based roughly on 7000 events for $`\sigma _x^{}`$ and $`\sigma _y^{}`$ 2000 events.

The expected beam parameters for CESR for the data discussed in this paper are given in Table 1.

All the parameters in Table 1 are given in the limit of zero bunch current. Our measurements were taken over a long time, about four months of CESR and CLEO running, and under many different machine conditions. Thus we expect only rough agreement with the parameters given in Table 1, but we should be sensitive to dynamic effects caused by the non-zero bunch currents. The horizontal $`\beta ^{}`$ and $`ϵ`$ are expected to depend primarily on the beam current, while the vertical $`ϵ`$ is expected to be influenced more by the measured luminosity over beam current (the specific luminosity). Therefore we present measurements in the horizontal direction as functions of bunch current, while measurements in the vertical are presented as functions of specific luminosity. We have previously observed that $`\beta _x^{}`$ is reduced by roughly a factor of two in colliding beam conditions due to the dynamic beta effect dbeta . Dynamic effects are also expected to produce a larger $`ϵ_y`$ and a smaller $`\beta _y^{}`$ and $`\beta _x^{}`$ than given in Table 1.

To be able to measure the angular spread of the colliding beams the momentum resolution must be well understood. We study CLEO’s momentum resolution with a simulated sample of $`e^+e^{}\mu ^+\mu ^{}`$ events with no momentum spread in the incoming beam particles. Recent examples show that, for high momentum tracks, the simulation of the tracking accurately models the data dstar hourglass . The resolution on the angle between the $`\mu ^+`$ and the $`\mu ^{}`$ depends on global properties of the event such as the pattern of hits in the transverse and longitudinal views of the SVX. More hits improves the resolution and the $`r\varphi `$ hits primarily effect the resolution transverse to the beam, while $`z`$ hits effect the longitudinal resolution. Thus the data are divided according to the SVX hit pattern. For each sub-sample, the resolution is parameterized as

$$A+B/Nhits=\mathrm{Resolution}$$
(5)
where $`Nhits`$ is the total number of tracking hits on both tracks. We determine the width of the momentum resolution as a function of $`Nhits`$, by selecting on a range of $`Nhits`$ and performing a one dimensional fit to a Gaussian plus a flat background to account for events not coming from beam collisions, mainly cosmic rays. To determine $`A`$ and $`B`$ for each SVX hit pattern, we fit the two dimensional distribution of resolution width versus $`Nhits`$. This gives us a set of $`A`$’s and $`B`$’s and combined with the observed SVX hit pattern and number of tracking hits we determine the resolution for each data event. The resolutions we extract on the angle between the $`\mu ^+`$ and $`\mu ^{}`$ is better than 1 mrad in both the vertical and horizontal directions for the highest number of tracking hits and SVX hits and up to 3 mrads for the worst case.

We then fit the two dimensional data distributions of angular width versus number of tracking hits for each SVX hit pattern with a width that is quadratic sum of the resolution piece determined from the simulation as described above and an extra contribution due to the underlying angular width of the colliding beams. We note that the means of the angular distributions are consistent with the known crossing angle, 2.0 mrads in the horizontal and none in the vertical, with an accuracy of better than a 0.01 mrads. The underlying width is then the weighted average of the fitted underlying widths over the SVX hit patterns. Systematic uncertainties on the extracted angular width of the colliding beams are applied by doing the data fits again with the values for the $`A`$ and $`B`$ parameters of the resolutions varied up or down in concert by one standard deviation based on the results from the simulation. Typically this represents a variation of 0.2–0.4 mrad. We also repeat the analysis using only the events with the largest number of tracking hits with the best resolution and see negligible changes.

We are not sensitive to the vertical width of the luminous region. Our resolution on the size of the luminous region is about $`30\mu `$m while our earlier observation showed that the underlying width is about $`7\mu `$m. We take the observed vertical width to be the resolution on the horizontal width, and this agrees well with the prediction of the simulation. For this reason, dynamic measurements of $`\beta _y^{}`$ cannot be extracted. A constant measured value, $`\beta _y^{}=(17910\pm 170)\mu `$hourglass , is used in this analysis to extract $`ϵ_y`$.

Figures 2

and 3 show the horizontal size and angular spread,

of the luminous region respectively, as functions of bunch current. These are after the resolutions discussed above have been taken out. Note that these are not the beam parameters but the result of the colliding beams. The size is a factor of $`\sqrt{2}`$ smaller than the beam size due to the overlap of the two beams producing the luminosity, and the angular size is a factor of $`\sqrt{2}`$ larger than the angular beam spread due to the incoherent collision of the two beams producing the luminosity. Figure 4 shows the width of the vertical angular

spread of the luminous region as a function of specific luminosity. From this data and algebraic manipulation of Equations 123 and 4 and the fixed $`\beta _y^{}`$ discussed above we extract $`\beta _x^{}`$, $`ϵ_x`$ and $`ϵ_y`$.

The extracted $`\beta _x^{}`$ is shown in Figure 5. This result agrees well with our

expectations based on a previous analysis done which observed the dynamic beta effect dbeta . The expectation was $`\beta _x^{}=1.2`$m reduced by roughly a factor of two. Here we see a $`\beta _x^{}0.6`$m. The extracted horizontal emittance as a function of specific luminosity is displayed in Figure 6. We expected an $`ϵ_x=0.21\mu `$m-rad

while we extracted an $`ϵ_x0.3\mu `$m. These results again agree well with our expectation. We do not clearly see any dynamic effects. We also expect an extracted vertical emittance of $`ϵ_y=0.0010\mu `$m-rad. We extract an $`ϵ_y0.0030\mu `$mrad. Our extracted value agrees with an earlier observation of $`ϵ_y=(0.0070\pm 0.0054\pm 0.0019)\mu `$m-rad hourglass  and with values calculated based on the observed luminosity and the best possible theoretical estimates of the other beam parameters. This extracted vertical emittance, shown in Figure 7

appears to have dynamic effects. When a fit to a first order polynomial is applied to this data we see a slope that is significant at only the 2.4 standard deviation level, which we do not claim as significant. We observe no visible dynamic effects on vertical emittance although the data are suggestive.

These methods can also be applied to observe the beam energy spread which results in a mismatch of the longitudinal momenta of the two tracks in these events. Our resolution is not precise enough to measure the energy spread, $`\mathrm{\Delta }E/E`$ is expected to be smaller than 0.0007 at CESR, resulting in an upper limit of 0.0015 at 90% C.L. for this data.

We have developed methods of measuring the emittances and betas of colliding $`e^+e^{}`$ beams in a non-destructive manner. The method depends on a very good understanding of the resolution of a charged particle tracking system observing the collision region. We demonstrate the methods using $`e^+e^{}\mu ^+\mu ^{}`$ events from the luminous region of CESR at the CLEO interaction point. We extract the beam parameters $`\beta _x^{}`$ and $`ϵ_x`$ as functions of bunch current and and $`ϵ_y`$ as a function of specific luminosity. The observations agree with our expectations and with measurements made by other methods. No dynamic effects are observed, although the data for $`ϵ_y`$ are suggestive of an increase in the emittance with specific luminosity.