text:
# Generating a random sink-free orientation in quadratic time

## 1. Introduction

A common problem is to select a random sample efficiently from a large collection of combinatorial objects. There are many reasons one may wish to do this. One is to obtain an approximate count: Jerrum and Sinclair \[JS\] showed that if one can generate nearly uniform samples, then for each $`\epsilon >0`$, one can obtain the cardinality of the collection to within a factor of $`1+\epsilon `$ with probability $`1\epsilon `$, in just a little more time. When counting the collection is #P-hard, as in the case of properly $`k`$-coloring a graph, this may be the only reasonable way to count, since it is unlikely that #P-hard counting problems can be solved exactly in polynomial time. Another reason to seek a sampling algorithm is that it may shed light on properties of the typical sample. For example, the analysis of typical spanning trees of Cayley graphs \[P, BLPS\] relies on two algorithms, the first developed by Aldous \[A\] and Broder \[B\] and the second by Wilson \[W\]. The analysis of phase boundaries in typical domino tilings of regions known as Aztec diamonds also relies on a sampling algorithm, known as domino shuffling \[CEP\]. Finally, sample generation may be a way of producing conjectures about the typical sample via simulation, when no theorem is known (for example, the results in \[CEP\] were initially discovered this way).

One common way to generate samples is Markov Chain Monte Carlo (MCMC). Here one finds an ergodic Markov chain whose equilibrium measure is the desired distribution $`\mu `$; then one runs the chain until the distribution is close to $`\mu `$. Constructing such a chain is usually easy (often when $`\mu `$ is uniform, there is a natural doubly stochastic transition matrix) and the hard part is knowing how long to run it. This may be established via eigenvalue bounds, or via coupling arguments or stopping times. In cases where the time bounds on the chain are established via coupling, it is often possible to improve on MCMC by using coupling from the past (CFTP) to obtain an exact sample rather than an approximate one \[PW1\].

In this note we consider the generation of a random sink-free orientation (SFO) of a finite undirected graph. Sink-free orientations were introduced by Bubley and Dyer \[BD\], who were motivated by an equivalence between counting them and counting satisfying assignments of Boolean formulas in conjunctive normal form in which each variable occurs at most twice (they call this problem Twice-SAT). Bubley and Dyer showed that counting sink-free orientations is #P-complete, so it is unlikely that an exact count can be obtained in polynomial time, and we must use approximate counting techniques based on nearly uniform sampling.

Bubley and Dyer give an MCMC algorithm that produces a sample whose distribution is within $`\epsilon `$ of uniform (in total variation) in time $`O(m^3\mathrm{log}(1/\epsilon ))`$, where $`m`$ is the number of edges. Huber \[H\] uses Bubley and Dyer’s analysis along with CFTP to produce an exact uniform sample in mean time $`O(m^4)`$. The purpose of this note is to improve the running time to $`O(nm)`$, where $`n`$ is the number of vertices. Instead of MCMC, we use a strong uniform time algorithm inspired by David Wilson’s cycle popping algorithm \[W\] for generating uniform directed spanning trees.

We now describe the problem and our results more precisely. Let $`G=(V,E)`$ be a finite undirected graph. We allow multiple edges and self-loops (but at most one self-loop per vertex, since multiple self-loops play no useful role in sink-free orientations). We define an $`n`$-cycle to be a ring of $`n`$ vertices $`v_0,\mathrm{},v_{n1}`$ with edges from $`v_i`$ to $`v_{i+1}`$ for each $`i`$ (taken modulo $`n`$; note that a $`1`$-cycle is a vertex with a self-loop), and an $`n`$-lollipop to be a path consisting of $`n`$ vertices and $`n1`$ edges, with a self-loop added at one end.

An orientation of an edge between vertices $`v`$ and $`w`$ is a mapping of the set $`\{\text{head},\text{tail}\}`$ onto $`\{v,w\}`$. Thus, a self-loop has only one orientation, but all other edges have two. To reverse the orientation of an edge, swap its head and tail. An orientation of $`G`$ is an orientation of each edge. A sink in an orientation is a vertex that is not the tail of any edge (a source is the opposite, i.e., not the head of any edge), and a sink-free orientation (SFO) of $`G`$ is an orientation that contains no sinks. If any connected component of $`G`$ is a tree, then $`G`$ has no SFO, and vice versa. Henceforth we restrict consideration to the class $`𝒮`$ of graphs in which no component is a tree. Let $`\mu _G`$ denote the probability measure assigning probability $`1/N`$ to each SFO of $`G`$, where $`N`$ is the total number of SFO’s of $`G`$.

Our algorithm, which we call “sink popping,” works as follows. Given a graph, orient the edges by independent, fair coin flips. If this orientation has no sinks, then it is the SFO we seek. Otherwise, choose any sink, and randomly re-orient each edge that points into the sink (i.e., all of its edges). We call this popping the sink, for reasons that will become clearer in the next section. Repeat until there are no more sinks.

We now state our main result.

###### Theorem 1.1.

For every graph $`G𝒮`$, sink popping terminates in finite time with probability $`1`$, regardless of how one chooses which sink to pop, and produces an output whose distribution is precisely $`\mu _G`$. The average number of sinks that must be popped is at most $`\left(\genfrac{}{}{0pt}{}{n}{2}\right)`$, where $`n`$ is the number of vertices of $`G`$, once again regardless of how one chooses which sink to pop. Equality holds only for the $`n`$-cycle or the $`n`$-lollipop. The expected number of times each particular vertex is popped is at most $`n1`$.

Sink popping is briefly mentioned at the end of \[PW2\], where the claims in the first sentence of the theorem are mentioned without detailed proof (and the running time is not analyzed).

To show that sink popping’s running time is $`O(nm)`$, we need to state the algorithm slightly more carefully. The subtle point is avoiding spending lots of time searching for sinks. We will keep a list of all sinks in the graph, and also a table showing the out-degree of every vertex. Generating these initially takes time proportional to the sum of the degrees of the vertices, or $`O(m)`$ time. Whenever we search for a sink to pop, we simply take the first sink from the list. When we pop the sink, we update the table to reflect the changes to its out-degree, and to those of its neighbors. It neighbors may have become sinks, in which case we append them to the list of sinks. (The purpose of the table is to let us easily see whether the neighbors have become sinks, without having to examine all their edges: in a complete graph, that would waste lots of time.) No sink can be annihilated except the one we popped, since no two sinks can share a common edge (it would have to point to both). Thus, each time we pop a sink at $`v`$, re-orienting its edges and updating the list and table requires time $`O(\mathrm{deg}(v))`$. By Theorem 1.1, the expected number of times $`v`$ is popped is $`O(n)`$, so the total expected number of operations is

$$O\left(\underset{v}{}n\mathrm{deg}(v)\right)=O(nm).$$

This time bound does not actually estimate the number of bit operations, but instead treats individual graph operations as units.

In the next section we give another description of the sink popping algorithm and explain its connection to cycle popping. We also state some further results about sink popping with arbitrary initial conditions. The third section contains proofs of the diamond and strong uniformity lemmas, which are analogous to the equivalent lemmas for cycle popping. The fourth section analyzes the running time. The fifth section derives some further facts about the running time. We conclude with some speculations and open questions.

## 2. Sink popping and cycle popping

Let $`H=(V,E)`$ be a finite, connected, directed graph, and let $`v`$ be any fixed vertex of $`H`$. A directed spanning tree of $`(H,v)`$ is a subset of edges so that every vertex other than $`v`$ has out-degree 1 and $`v`$ has out-degree 0. Wilson \[W, PW2\] invented the following algorithm, known as “cycle popping,” for generating a uniform random directed spanning tree. For each $`wV\{v\}`$ and $`k0`$, let $`X_{w,k}`$ be a random edge leading out of $`w`$, chosen uniformly from among all edges leading out of $`w`$. Let these be independent as $`w`$ and $`k`$ vary. For fixed $`w`$, imagine the collection $`\{X_{w,k}:k0\}`$ as a stack with $`X_{w,0}`$ on top. Initially, look at the collection $`\{X_{w,0}:wV\{v\}\}`$, that is, consider the collection $`\{X_{w,f(w)}\}`$ with $`f0`$. If these form a directed spanning tree, stop and set the sample equal to it. If not, there must be a cycle in this collection. Choose a cycle (it doesn’t matter which), and increment $`f(w)`$ by 1 for each $`w`$ in the cycle. (Imagine popping these edges off the stack so the next element of each stack is now on top.) If the collection $`\{X_{w,f(w)}:wV\{v\}\}`$ is now a directed spanning tree, stop and return this for your sample, otherwise continue popping until you do stop. Wilson showed that the set of cycles popped does not depend on which you choose to pop when you have a choice, and that the algorithm stops almost surely at a directed spanning tree with uniform distribution.

We can describe sink popping in similar terms, which will be useful in the proof of Theorem 1.1. Let $`G=(V,E)`$ be a finite undirected graph in $`𝒮`$ and let $`\mathrm{\Omega }=\mathrm{\Omega }_0^{}`$, where $`\mathrm{\Omega }_0`$ is the set of orientations of $`G`$, i.e., $`\mathrm{\Omega }`$ consists of sequences of orientations of $`G`$. We endow $`\mathrm{\Omega }`$ with the $`\sigma `$-field $``$ generated by the coordinate functions $`X_{e,k}`$ for $`eE(G)`$ and $`k0`$, which specify the orientation of $`e`$ in the $`k`$-th orientation in the sequence. Endow $`\mathrm{\Omega }`$ with the probability measure $`𝐏`$ under which the coordinate functions are independent and each equally likely to yield either orientation. The intuition is that $`\{X_{e,k}:k0\}`$ represents a stack of arrows under the edge $`e`$. Define a random function $`f:E\times `$ as follows. Let $`f(e,0)=0`$ for all $`e`$. Given $`f(e,k)`$ for all $`e`$, define $`f(,k+1)`$ inductively: If the collection $`\{X_{e,f(e,k)}:eE\}`$ is an SFO, then set $`f(e,k+1)=f(e,k)`$ for all $`e`$. If not, choose a sink $`v_kV`$ arbitrarily, i.e., a vertex $`v_k`$ for which all edges $`e`$ incident to it are oriented toward it by the orientation $`X_{e,f(e,k)}`$. Let $`f(e,k+1)=f(e,k)`$ for $`e`$ not incident to $`v_k`$, and $`f(e,k+1)=f(e,k)+1`$ for $`e`$ incident to $`v_k`$. The dependence of $`f`$ on the choice rule (for choosing $`v_k`$, if there are several sinks) is suppressed in the notation, as is the dependence on the choice of $`\omega \mathrm{\Omega }`$ via the variables $`\{X_{e,k}\}`$. Intuitively, $`f(e,k)`$ is the original depth of the arrow under $`e`$ now at the top of the stack at time $`k`$. Say that $`v_k`$ is the sink popped at time $`k`$, and let $`\tau =\mathrm{min}\{k:\{X_{e,f(e,k)}\}\text{ is an SFO}\}`$ be the number of pops before an SFO is obtained (conceivably $`\tau =\mathrm{}`$), and $`\eta =\eta (\omega ,\text{choice rule})`$ denote the resulting SFO (if any).

Except for its last sentence, Theorem 1.1 is established by showing that $`\tau <\mathrm{}`$ with probability $`1`$, the law of $`\{X_{e,f(e,\tau )}\}`$ is precisely $`\mu _G`$, and $`𝐄\tau \left(\genfrac{}{}{0pt}{}{n}{2}\right)`$, with equality in and only in the cases indicated. The first two lemmas are analogous to those used by Wilson \[W\] in establishing the validity of the cycle popping algorithm, and the third is the running time analysis.

###### Lemma 2.2 (Diamond lemma).

The number of pops $`\tau \mathrm{}`$ in a maximal popping sequence is independent of the choice rule, as is the multiset $`\{v_k:0k<\tau \}`$. If $`\tau <\mathrm{}`$ then the resulting SFO $`\eta `$ is also independent of the choice rule.

The name “diamond” is meant to remind the reader that moving from the top of a diamond to the bottom by going southeast then southwest is equivalent to going southwest then southeast. This terminology comes from the article \[E\].

###### Lemma 2.3 (Strong uniform time).

Let $`N`$ be the number of SFO’s of $`G`$. Then for each $`k0`$, and each SFO $`\eta `$,

$$𝐏(\tau =k,\{X_{e,f(e,\tau )}\}=\eta )=\frac{𝐏(\tau =k)}{N}.$$

In other words, $`\tau `$ is a strong uniform time.

###### Lemma 2.4.

If $`G𝒮`$ has $`n`$ vertices, then $`𝐄\tau \left(\genfrac{}{}{0pt}{}{n}{2}\right)`$, with equality only for the $`n`$-cycle and the $`n`$-lollipop.

We conclude this section by stating two results that shed further light on the running time of the popping algorithm.

###### Proposition 2.5.

The distribution of $`\tau `$ for the $`n`$-cycle is exactly the same as the distribution for the $`n`$-lollipop.

###### Proposition 2.6.

Let $`G𝒮`$ be any graph with $`n`$ vertices and let $`_0`$ be the $`\sigma `$-field generated by the variables $`\{X_{e,0}\}`$. Then the conditional mean running time $`𝐄(\tau |_0)`$ is always bounded by $`n(n1)`$, and the only case to achieve this is an $`n`$-lollipop with all edges oriented opposite to their orientation in the unique SFO.

## 3. Strong uniformity

We first establish deterministic facts holding for every sample $`\omega \mathrm{\Omega }`$. Say that a sequence $`v_0,\mathrm{},v_{k1}`$ with $`k\mathrm{}`$ is a maximal popping sequence for $`\omega `$ if it is legal (i.e., only sinks are popped) and cannot be extended to larger $`k`$ (thus if $`k<\mathrm{}`$ it results in an SFO). Note that if $`k=\mathrm{}`$, we do not mean our notation to suggest that $`v_0,v_1,\mathrm{}`$ is followed by a final term $`v_\mathrm{}1`$; instead, $`v_0,\mathrm{},v_\mathrm{}1`$ denotes the infinite sequence $`v_0,v_1,\mathrm{}`$, with no final term.

Let $`f(e,k)`$ denote the function $`f(e,k,\omega ,𝐯)`$ where $`\omega `$ is a sample point and $`𝐯`$ is a specified legal sequence of pops of length at least $`k`$. Define an equivalence relation on finite sequences $`v_0,\mathrm{},v_{k1}`$ of vertices of $`G`$ by calling two sequences equivalent if one can be changed to the other by a sequence of transpositions of pairs of vertices $`(v_i,v_{i+1})`$ that are not neighbors in $`G`$. (Note that such a transposition does not change whether a sequence is a legal popping sequence.) The following lemma is useful, though obvious.

###### Lemma 3.7 (Deterministic strong Markov property).

Given an integer $`j`$ and vertices $`v_0,\mathrm{},v_{j1}`$, let $`\omega `$ be any initial configuration for which $`v_0,\mathrm{},v_{j1}`$ is a legal popping sequence and let $`\omega ^{}`$ and $`\omega `$ be related by

$$X_{e,k}(\omega ^{})=X_{e,k+f(e,j)}(\omega ).$$

That is, $`\omega ^{}`$ looks like $`\omega `$ after $`v_0,\mathrm{},v_{j1}`$ are popped. Then the following deterministic strong Markov property (DSMP) holds. For any $`k\mathrm{}`$, the set of sequences $`\{v_{j+i}:0ik1\}`$ for which $`v_0,\mathrm{},v_{j+k1}`$ is a legal popping sequence for $`\omega `$ is the same as the set of legal popping sequences of length $`k`$ for $`\omega ^{}`$. If $`v_0,\mathrm{},v_{j+k1}`$ is maximal for $`\omega `$ then $`v_j,\mathrm{},v_{j+k1}`$ is maximal for $`\omega ^{}`$, and leaves the same SFO (if $`k<\mathrm{}`$).

Extend the definition of equivalence to infinite sequences by saying that $`v_0,v_1,\mathrm{}`$ is equivalent to $`w_0,w_1,\mathrm{}`$ if, by a sequence of transpositions as above applied to $`v_0,v_1,\mathrm{}`$, one can transform $`v_0,v_1,\mathrm{}`$ so that arbitrarily long initial segments of it match those of $`w_0,w_1,\mathrm{}`$. In particular, this implies that the multisets $`\{v_k\}`$ and $`\{w_k\}`$ are the same.

Let $`l(\omega )`$ denote the minimal length of a maximal popping sequence for $`\omega `$.

###### Lemma 3.8.

The set of maximal popping sequences is an equivalence class.

###### Proof.

Let $`v_0,\mathrm{},v_{k1}`$ be a legal popping sequence for $`\omega `$ with $`k<\mathrm{}`$, and let $`w_0,\mathrm{},w_{k1}`$ be obtained from $`v_0,\mathrm{},v_{k1}`$ by transposing $`v_i`$ and $`v_{i+1}`$ which are not neighbors in $`G`$. Suppose $`i=0`$. Since the edges incident to $`v_0`$ are disjoint from the edges incident to $`v_1`$, we may apply the DSMP to $`v_0,v_1`$ and to $`v_1,v_0`$ and see that $`w_0,\mathrm{},w_{k1}`$ is legal as well and maximal if $`v_0,\mathrm{},v_{k1}`$ is. If $`i>0`$, first apply the DSMP to $`v_0,\mathrm{},v_{i1}`$ and then use the same argument. This shows that equivalent sequences are either both maximal popping sequences or neither. (The case of $`k=\mathrm{}`$ is trivial, since infinite popping sequences are automatically maximal.)

To prove the lemma, we induct on $`l(\omega )`$, and then deal with the case of $`l(\omega )=\mathrm{}`$. It is clear when $`l=0`$. Assuming the lemma for $`l(\omega )<L`$, let $`l(\omega )=L`$ with maximal popping sequence $`v_0,\mathrm{},v_{L1}`$. Let $`w_0,\mathrm{},w_{k1}`$ be any other maximal popping sequence. If $`w_0=v_0`$, applying the DSMP and the induction hypothesis completes the induction. If not, then consider the least $`i`$ for which $`v_i=w_0`$, if any. When we pop $`v_0,\mathrm{},v_{L1}`$, the orientation $`\{X_{e,f(e,j)}:eE\}`$ has a sink at $`w_0`$ for each $`j<i`$, since the sink at $`w_0`$ exists until one of its edges is popped and no other sink can contain any such edge until $`w_0`$ is popped. Thus $`i`$ exists and $`v_j`$ cannot be a neighbor of $`v_i`$ for $`j<i`$. Hence, we can move $`v_i`$ to the first position by a sequence of adjacent transpositions with non-neighbors in $`G`$. We have seen that the resulting sequence $`v_i,v_0,v_1,v_2,\mathrm{},v_{i1},v_{i+1},\mathrm{},v_{L1}`$ is a maximal popping sequence. Now apply the DSMP and the induction hypothesis to conclude that $`w_0,\mathrm{},w_{k1}`$ is equivalent to $`v_i,v_1,\mathrm{}v_{i1},v_{i+1},\mathrm{},v_{L1}`$, and thus to $`v_0,\mathrm{},v_{L1}`$.

All that remains is the case of $`l(\omega )=\mathrm{}`$, i.e., the case when all maximal popping sequences are infinite. Given two such sequences $`v_0,v_1,\mathrm{}`$ and $`w_0,w_1,\mathrm{}`$, the argument from the previous paragraph shows that we can transform $`v_0,v_1,\mathrm{}`$ so that its first element is $`w_0`$. Now applying the DSMP shows that we can bring arbitrarily long initial segments into agreement, which is the definition of equivalence. ∎

###### Proof of the diamond lemma.

From the previous lemma, we know that $`\tau =l`$, so $`\tau `$ is independent of the choice rule. Furthermore, since all maximal popping sequences are equivalent, the multisets of popped vertices are the same. The assertion about SFO’s follows because the SFO depends only on which vertices were popped. ∎

###### Proof of strong uniform time.

We prove by induction that for any SFO $`\eta `$ and any finite sequence $`v_0,\mathrm{},v_{k1}`$, the following event has probability $`2^{m_{i=0}^{k1}\mathrm{deg}_0(v_i)}`$, where $`\mathrm{deg}_0`$ means the degree not counting self-loops: $`\tau =k`$, and $`v_0,\mathrm{},v_{k1}`$ is a legal popping sequence for $`\omega `$, and $`\eta (\omega )=\eta _0`$. This is vacuously true when $`k=0`$. Now the probability that the singleton $`v_0`$ is a legal pop is $`2^{\mathrm{deg}_0(v_0)}`$, so applying the DSMP we see that the probability of a maximal popping sequence $`v_0,v_1,\mathrm{},v_{k1}`$ with $`\eta =\eta _0`$ is

$$2^{\mathrm{deg}_0(v_0)}2^{m_{i=1}^{k1}\mathrm{deg}_0(v_i)},$$

which completes the induction.

To find the probability of both $`\tau =k`$ and $`\eta =\eta _0`$ (with no restrictions on the popping sequence), we must sum this probability over all equivalence classes of potential popping sequences of length $`k`$. We sum over equivalence classes to avoid double counting, since for any given $`\omega `$, Lemma 3.8 tells us that the set of maximal popping sequences is an equivalence class. Since neither the summand nor the set of sequences depends on $`\eta _0`$, we have proved the lemma. ∎

## 4. Analysis of the running time

We still have not shown that $`\tau `$ is almost surely finite. While this may appear obvious from some kind of Markov property, the choice rule makes things sticky and we find it easiest to conclude this from the existence of a finite upper bound on the expected run time. To bound $`\tau `$ we make repeated use of the following monotonicity principle. We let $`Q(G,v)`$ denote the random number of times $`v`$ is popped in a maximal popping sequence (possibly $`\mathrm{}`$), which, by the diamond lemma, is well defined.

###### Lemma 4.9 (Monotonicity).

Fix $`G𝒮`$ and let $`H𝒮`$ be a subgraph of $`G`$, that is, $`V(H)V(G)`$ and $`E(H)E(G)`$. For $`vV(H)`$,

$$𝐄Q(H,v)𝐄Q(G,v).$$

###### Proof.

This is proved by stochastic domination: we run sink popping simultaneously on $`H`$ and $`G`$, using the same stacks for edges common to both graphs. Every legal popping sequence on $`G`$ restricts to a legal popping sequence on $`H`$ as well, so under this coupling $`Q(H,v)Q(G,v)`$ always. ∎

Remark: Additionally, we see that equality occurs only when no SFO on $`G`$ can require further popping of $`v`$ on $`H`$.

###### Proposition 4.10.

If $`G`$ is an $`n`$-cycle, then $`𝐄\tau =\left(\genfrac{}{}{0pt}{}{n}{2}\right)`$. Furthermore, conditioned on starting with $`j`$ edges oriented clockwise and $`nj`$ counterclockwise, the expected value of $`\tau `$ is $`2j(nj)`$.

###### Proof.

At any time, some of the arrows point clockwise and others counterclockwise. Let $`Y_k`$ be the number of arrows pointing clockwise at time $`k`$. Popping at any vertex causes two opposite pointing arrows to be replaced by two random arrows. Thus $`Y_{k+1}`$ has the distribution of $`Y_k+Z`$ where $`𝐏(Z=1)=𝐏(Z=1)=1/4`$ and $`𝐏(Z=0)=1/2`$. Therefore $`\{Y_k:k0\}`$ is a simple random walk with delay probability of $`1/2`$ absorbed at $`0`$ and $`n`$. The expected absorption time from $`j`$ is twice that for simple random walk, and thus is $`2j(nj)`$; see equation (3.5) on page 349 of Feller \[F\]. Hence $`𝐄\tau =2𝐄Y_0(nY_0)`$, which is twice the expected number of ordered pairs of edges where the first is initially clockwise and the second initially counterclockwise. There are $`n(n1)`$ ordered pairs of distinct edges, each having these orientations with probability $`1/4`$, so $`𝐄\tau =n(n1)/2`$. ∎

###### Corollary 4.11.

Let $`𝒮_0`$ denote the class of graphs in which every vertex is in some cycle. For $`G𝒮_0`$ and $`vV(G)`$,

$$𝐄Q(G,v)(n1)/2.$$

Equality holds for all $`v`$ if and only if $`G`$ is an $`n`$-cycle.

###### Proof.

Fix $`G`$ and $`v`$ and let $`H`$ be a cycle containing $`v`$. By monotonicity, $`𝐄Q(G,v)𝐄Q(H,v)`$ which is at most $`(n1)/2`$ by Proposition 4.10 and symmetry. In fact, $`𝐄Q(H,v)`$ is strictly less than $`(n1)/2`$ unless $`H`$ is an $`n`$-cycle. By the remark following the proof of the monotonicity lemma, the inequality $`𝐄Q(G,v)𝐄Q(H,v)`$ is strict unless no SFO on $`G`$ can require further popping of $`v`$ on $`H`$. In our case, $`H`$ is an $`n`$-cycle, and $`G`$ is an $`n`$-cycle with some chords or self-loops added. Then (assuming $`G`$ is not an $`n`$-cycle), there is always an SFO on $`G`$ that does not restrict to an SFO on $`H`$: if $`G`$ has a self-loop, one can choose an SFO on $`G`$ such that $`H`$ has a sink there; if $`G`$ has a chord, one can use the chord to create a short circuit across $`H`$ giving a cycle of length less than $`n`$, orient this cycle in a loop, and orient the other edges in $`G`$ towards the cycle. If $`v`$ is a sink in the restriction to $`H`$ of such an SFO on $`G`$, then strict inequality holds for $`v`$ (because with positive probability, sink popping on $`G`$ will produce this SFO, and $`v`$ will still need to be popped in $`H`$). ∎

###### Lemma 4.12.

For every $`G𝒮`$ with $`n`$ vertices, and each $`vV(G)`$,

$$𝐄Q(G,v)n1,$$

and equality holds only when $`v`$ is the vertex furthest from the self-loop in an $`n`$-lollipop.

###### Proof.

We induct on $`G`$. The base step is $`G𝒮_0`$, which is immediate from the previous corollary. Assume for induction that the conclusion holds for all subgraphs of $`G`$. There are three cases other than the base step.

Case 1. $`G`$ is not connected. Then the result follows from the induction hypothesis and the monotonicity lemma applied to the component $`H`$ of $`G`$ containing $`v`$. Equality never occurs.

If $`G`$ is connected and not in $`𝒮_0`$, then $`G`$ must contain an isthmus, i.e., an edge whose removal disconnects $`G`$.

Case 2. Some edge $`e`$ disconnects $`G`$ into two components both in $`𝒮`$. Again the result follows from the induction hypothesis applied to the component $`H`$ of $`G\{e\}`$ containing $`v`$, and equality never occurs.

Case 3. $`G`$ has an isthmus and removal of any isthmus always leaves a component that is a tree. Then $`G`$ has a leaf $`z`$. If $`vz`$ then the result follows immediately from monotonicity with $`H=G\{z\}`$. If $`v`$ is the only leaf, then let $`w`$ be its neighbor. Choose a popping order that pops $`v`$ whenever possible, and otherwise executes any choice rule for sink popping on $`H:=G\{v\}`$. Initially there is a $`1/2`$ chance that $`v`$ is a sink, in which case it is popped a mean $`2`$ geometric number of times until the edge $`vw`$ points to $`w`$. Then, each time $`w`$ is popped, the probability is $`1/2`$ that this edge is reversed, in which case it takes another mean $`2`$ geometric number of pops to reverse it again. Thus

$$𝐄Q(G,v)=1+𝐄Q(H,w).$$

By induction, this is at most $`1+(n2)`$. Equality occurs for $`v`$ in $`G`$ if and only if it occurs for $`w`$ in $`H`$, so we see by induction that it holds only at the end of a lollipop. ∎

###### Proof of Lemma 2.4.

We prove the lemma by induction, following the pattern of the last proof. The base step is $`G𝒮_0`$, in which case the lemma follows from Corollary 4.11. In the cases 1 and 2 of the induction, if $`G`$ is disconnected or the union of two graphs in $`𝒮`$ along an added edge, the result is again immediate from the subadditivity of the function $`n\left(\genfrac{}{}{0pt}{}{n}{2}\right)`$ and monotonicity. Finally, if $`G`$ has a leaf $`v`$, we set $`H:=G\{v\}`$ and observe that the number of pops $`\tau _G`$ and $`\tau _H`$ on $`G`$ and $`H`$ respectively are related by $`\tau _G=\tau _H+Q(G,v)`$. Thus

$$𝐄\tau _G=𝐄\tau _H+𝐄Q(G,v)\left(\genfrac{}{}{0pt}{}{n1}{2}\right)+(n1)=\left(\genfrac{}{}{0pt}{}{n}{2}\right).$$

By the previous lemma, the last inequality is strict unless $`H`$ is an $`n`$-lollipop and its vertex of degree 1 is the neighbor of $`v`$ in $`G`$. This completes the induction. ∎

###### Proof of Theorem 1.1.

The theorem follows immediately from combining Lemmas 2.2, 2.3, 2.4, and 4.12. ∎

## 5. Further proofs

The $`n`$-cycle and $`n`$-lollipop have the worst mean run times. Here we prove Proposition 2.5, namely that the run time distributions are in fact identical.

###### Proof of Proposition 2.5.

Number the vertices of the $`n`$-lollipop $`0,\mathrm{},n1`$ with $`0`$ being the leaf. Always pop the sink with lowest number. Let $`Y_k`$ denote the sink popped at time $`k`$. Clearly $`Y_0`$ is $`1`$ plus a mean $`2`$ geometric random variable, with the proviso that a value of $`n1`$ or higher represents the terminal state in which no sink needs to be popped. Let $`_k`$ be the $`\sigma `$-field generated by $`Y_0,\mathrm{},Y_{k1}`$. We claim that $`\{Y_k:k0\}`$ is a time-homogeneous Markov chain with respect to $`\{_k\}`$ and that from any state $`j>0`$ its increments are $`2`$ plus a mean $`2`$ geometric, jumping to the terminal state if it reaches $`n1`$ or greater, and from state $`0`$ the same thing with $`2`$ replaced by $`1`$ (thus the jump from $`0`$ is resampled if it hits $`1`$). All that is needed to check this is an inductive verification that the orientations of edges between vertices of higher index than $`Y_k`$ are conditionally i.i.d. fair coin flips given $`_k`$, which is straightforward.

Now we show that the running time on an $`n`$-cycle is also equal to the time for a random walk to hit at least $`n1`$ when its increments are $`2`$ plus a mean $`2`$ geometric, resampled if it hits $`1`$. At time $`k`$, let $`Y_k`$ denote the least index of a sink when the edge from $`n1`$ to $`0`$ is oriented toward $`0`$ and $`n1`$ minus the greatest index of a sink when the edge is oriented toward $`n1`$. In other words, this quantity is the distance from the head of the $`0,n1`$ edge to the nearest sink in that direction. We always choose the pop that sink. The only time the $`0,n1`$ edge can change orientations is when $`Y_k=0`$, in which case $`Y_{k+1}`$ will be $`1`$ plus a mean $`2`$ geometric; when $`Y_k>0`$ verification of the conditional increment is trivial. The stopping rule is, again, that one must jump to $`n1`$ or greater, and $`Y_0`$ has the right distribution for the same reason as before, so the sequence has the same distribution. ∎

Our final result deals with the run time started from an arbitrary state, that is, the conditional distribution of $`\tau `$ given $`_0`$ (as defined in Proposition 2.6). While this quantity is a hidden variable as far as users of the algorithm are concerned, it has relevance to the distribution of the run time, as well as having some intrinsic interest. We begin again with a result on the $`n`$-cycle.

###### Proposition 5.13.

Let $`G`$ be an $`n`$-cycle. Then for every $`vV(G)`$,

$$𝐄(Q(G,v)|_0)3n/4,$$

with equality if and only if $`n`$ is even and all edges are oriented along the direction of shortest travel to $`v`$.

###### Proof.

Number the vertices $`0,\mathrm{},n1`$ mod $`n`$. We first establish that the discrete Laplacian of $`𝐄(Q(G,v)|_0)`$ depends on the initial orientation of $`G`$ via

$$𝐄\left[Q(G,v)\frac{Q(G,v+1)+Q(G,v1)}{2}|_0\right]=\{\begin{array}{cc}1\hfill & \text{if }v\text{ is a sink,}\hfill \\ 1\hfill & \text{if }v\text{ is a source, and}\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$
(5.1)
To see this, choose any popping order and let $`Y(v,k)`$ denote the in-degree of $`v`$ at time $`k`$, that is, the number of $`eE(G)`$ adjacent to $`v`$ for which $`X_{e,f(e,k)}`$ is oriented toward $`v`$. Then, conditionally on anything up to time $`k`$,

$$𝐄Y(v,k+1)=𝐄Y(v,k)𝐏(v_k=v)+\frac{𝐏(v_k=v+1)+𝐏(v_k=v1)}{2},$$

since any pop at $`v`$ reduces the expected in-degree by mean 1 and any pop at a neighbor of $`v`$ increases it by mean $`1/2`$. Summing over $`k`$, conditioning on $`_0`$ and using $`Y(v,\tau )1`$ proves (5.1).

From Proposition 4.10 we know that $`𝐄(\tau |_0)=2Y_0(nY_0)`$ where $`Y_0`$ is the number of initial clockwise arrows (edges oriented from $`i+1`$ to $`i`$ mod $`n`$ for some $`i`$). This, along with (5.1), determines $`𝐄(Q(G,)|_0)`$, since the difference of any two candidates for this function would be a harmonic function on the cycle, and hence constant.

In general,

$$𝐄(Q(G,v)|_0)=\frac{k}{n}\left(1+3n2k\frac{2}{k}\underset{j=1}{\overset{k}{}}a_j\right),$$
(5.2)
if there are $`k`$ clockwise edges pointing from $`v+a_j`$ to $`v+a_j1`$ for a set $`\{a_1,\mathrm{},a_k\}\{1,\mathrm{}n\}`$ (addition taken mod $`n`$). To prove this formula, we need only prove that the right hand side satisfies the two properties that characterize the left hand side. The sum over all $`v`$ (i.e., $`𝐄(\tau |_0)`$) is easy, since it equals

$$k(1+3n2k)\frac{2}{n}\underset{j=1}{\overset{k}{}}\underset{i=1}{\overset{n}{}}i,$$

which does indeed simplify to $`2k(nk)`$. To check that the right hand side of (5.2) works in (5.1), we proceed as follows. Let $`f(v)`$ be the right hand side of (5.2). Then

$$g(v)\stackrel{\text{def}}{=}f(v+1)f(v)=\frac{2}{n}\left(k+n\delta _{\mathrm{min}\{a_j\},1}\right),$$

where $`\delta `$ is the Kronecker delta. Hence,

$$\frac{g(v)g(v1)}{2}=\{\begin{array}{cc}1\hfill & \text{if }v\text{ is a sink,}\hfill \\ 1\hfill & \text{if }v\text{ is a source, and}\hfill \\ 0\hfill & \text{otherwise,}\hfill \end{array}$$

as desired.

Equation (5.2) makes it easy to see when $`𝐄(Q(G,v)|_0)`$ is maximized: that can occur only when $`\{a_1,\mathrm{},a_k\}=\{1,\mathrm{},k\}`$, in which case $`𝐄(Q(G,v)|_0)`$ equals $`3n(k/n)(1k/n)`$. This quantity is bounded above by $`3n/4`$, with equality if and only if $`n=2k`$. ∎

###### Proof of Proposition 2.6.

Induct again, as in the proof of Lemma 4.12 and the main theorem. Simultaneously, we show by induction that $`𝐄(Q(G,v)|_0)2(n1)`$, with equality only for the leaf of an $`n`$-lollipop and initial conditions $`X_{e,0}`$ all pointing toward $`v`$. Note that the previous proposition proves this bound for all $`G𝒮_0`$, if we use monotonicity (which also holds conditioned on the initial orientation).

For the base case, $`G𝒮_0`$ and Proposition 5.13 shows that in fact $`𝐄(\tau |_0)<3n^2/4`$. There is strict inequality because equality in Proposition 5.13 cannot hold simultaneously for all vertices in a cycle. It follows that $`𝐄(\tau |_0)<n(n1)`$ unless $`n3`$. The cases with $`n3`$ are easily dealt with: those with $`n=2`$ are trivial, and for $`n=3`$ the worst case contains a $`3`$-cycle, which can be analyzed using the sharper bounds in the proof of Proposition 5.13.

When $`G`$ is not connected or is the union of two graphs in $`𝒮`$ along an added isthmus, the $`n(n1)`$ bound is immediate from subadditivity of $`n(n1)`$ and monotonicity, and the $`2(n1)`$ bound follows from monotonicity. Finally, when $`G`$ has a leaf $`v`$, set $`H:=G\{v\}`$ as before. This time, in the worst case we know that $`v`$ is a sink initially, so $`𝐄(Q(G,v)|_0)`$ is bounded by $`2+𝐄(Q(H,w)|_0)`$. This verifies the conclusion that $`𝐄(Q(G,v)|_0)2(n1)`$, and adding this to $`𝐄\tau _H`$ gives, by induction, at most $`(n1)(n2)+2(n1)=n(n1)`$, which completes the proof of the upper bound; the conditions for equality are clear from the proof. ∎

## 6. Questions

It is tempting to view both cycle popping and sink popping as special cases of what might be called “partial rejection sampling:” to generate a random structure, choose a random candidate, and if it has any flaws, locally rerandomize until it is flawless. Does partial rejection sampling apply to other natural combinatorial problems? Can one develop a general theory? Note that Fill and Huber’s randomness recycler \[FH\] also uses the idea of rejecting only part of a structure, although in a different way.

Cycle popping was applied to the study of random spanning trees on $`^d`$, as well as some more general graphs, in \[BLPS\]. It would be interesting if sink popping could be used similarly. Do random sink-free orientations on $`^d`$ exhibit any interesting or surprising structure?

## Acknowledgements

We thank David Wilson and Peter Winkler for helpful discussions, and the anonymous referee for useful comments on the manuscript.


format: text

text:
# The dark matter crisis

## Introduction

Most dark matter candidates cannot be distinguished by observations on large scales. Although the observed universe appears consistent with hierarchical structure formation, this still leaves a wide range of potential dark matter candidates. For example, the clustering properties of galaxies, abundances of rich clusters or even halo masses and sizes are all very similar in universes with matter density dominated by cold dark matter, warm dark matter or collisional dark matter. We therefore seek tests of the nature of the dark matter that are sensitive to its interaction properties and small scale power which manifests itself on non-linear scales.

Cold dark matter (CDM) halos form via a complicated sequence of hierarchical mergers that lead to a global structure set primarily by violent relaxation. Numerical simulations have played an important role in determining the shape and scaling of CDM halo profiles that have subsequently lead to new observational tests of the model. A single functional form can fit CDM halos from a mass scale of $`10^7M_{}`$$`10^{15}M_{}`$, where the density at a fixed fraction of the virial radius is higher in lower mass halos and the central profiles have steep singular cusps (c.f. Dubinski & Carlberg 1991, Warren etal 1992, Navarro etal 1996, Fukushige & Makino 1997, Moore etal 1998, Jing 2000 etc).

Galaxy clusters form via a similar process as individual galaxies, however most of the galactic fragments that formed clusters have survived the hierarchical growth, whereas on galaxy scales we find little trace of the merging hierarchy. Only a dozen satellites orbit the Milky Way, whereas a thousand satellites orbit within the Coma cluster. Numerical simulations of CDM halo formation have revealed that the abundance of dark matter subhalos within a galaxy is the same as found within a scaled galaxy cluster (Moore etal 1999, Klypin etal 1999). There are two solutions to this problem: (i) CDM is incorrect and the nature of the dark matter suppresses the formation of substructure halos, (ii) CDM is correct and the dark matter satellites of the Milky Way are present but only a few percent of them formed stars. Here I focus on the latter possibility.

## I Dark matter substructure

Figure 1. The distribution of dark matter with a CDM “Local Group”candidate. This is a binary pair of dark matter halos at a redshift z=0, separated by 1 Mpc and infalling at 100km/s. The large halos have virial masses of $`2\times 10^{12}M_{}`$, and with a particle mass of $`10^6M_{}`$ they are resolved with over $`10^6`$ particles and 0.5 kpc force resolution. The grey scale represents the local density of dark matter – there are over 2000 dark matter satellites with circular velocity larger than 10 km/s.

Figure 2. The cumulative number of satellites within the virial radius of the Milky Way (open circles) and within the two CDM halos from Figure 1 (open triangles). Here I have taken $`v_{peak}=210`$ km/s for the Milky Way, however the CDM contribution to the Milky Way can be constrained to lie in the range 130 – 160 km/s once the baryonic component has been considered (Moore etal 2001). The dotted curves show the effects of this correction. The arrows show a correction for converting central velocity dispersion to $`v_{sat}`$. The thick solid curve shows the distribution of CDM satellites that could form stars before the universe is re-ionised.

Interpreting the observations

Figure 2 shows the cumulative distribution of subhalos within the high resolution Local Group halos of Figure 1. The open circles show the observed distribution for the Milky Way satellites, where I have normalised the distribution using $`v_{peak}=210`$ km/s. However, baryons dominate the central region of the Galaxy and subtracting the contribution from the disk and bulge gives the maximum allowed CDM halo that has $`v_{peak}=160`$ km/s (Moore etal 2001). A minimum value of $`v_{peak}=130`$ km/s is required for the Galactic CDM halo to be massive enough to cool the observed mass of baryons. Figure 2 shows the effect of this correction.

Simon White has pointed out that the velocity dispersion of the dSph’s are measured well within the cores of their dark matter halos (White 2000). We originally assumed isotropic orbits and isothermal potentials to derive $`v_{sat}`$ from observations of the 1d central velocity dispersion (Moore etal 1999). CDM halos have central density profiles flatter than $`r^2`$ therefore one expects the velocity dispersion to drop in the inner region. M87 provides a good example of this. This galaxy lies at the center of the Virgo cluster and has a central velocity dispersion of $`350`$ km/s whereas the cluster has a global value that is a factor of two larger. If we assume that the dSph’s are similar to M87, then the correction should scale roughly as the concentration parameter. Since $`c_{M87}/c_{dSph}0.40.5`$, then we expect the maximum correction to $`v_{peak}`$ to be an increase of 50% over our quoted values. This correction is indicated by the arrows in Figure 2 which brings the observed data into good agreement with a crude model for re-ionisation discussed later.

## II Feedback

Figure 3. The histogram shows the cumulative distribution of absolute magnitudes of the 11 Galactic satellites. The curves show predictions from the Durham semi-analytic models of galaxy formation (Cole etal 2000) varying the parameter $`\alpha `$ that controls the efficiency of feedback. This can be tuned to give the correct distribution of satellite luminosities which in this case would lie somewhere in between the blue and pink curves.

Feedback is an essential component of galaxy formation within CDM models. It is invoked primarily to flatten the luminosity function given the steep mass function of CDM halos. For example, the faint end of the luminosity function in the Local Group is flat over a range that is about 10 magnitudes fainter than $`M_{}`$. By varying the efficiency of feedback with halo mass, it is possible to get a reasonably flat luminosity function as Figure 3 demonstrates. The parameter $`\alpha `$ controls how much gas is ejected from dark matter halos of a given circular velocity and allows one to form systematically less stars in smaller mass halos.

The problem with a uniform feedback scheme is that the mass to light ratios of galaxies will increase rapidly for fainter galaxies. Thus we find that a satellite halo with absolute magnitude $`M_B=10`$ is predicted to have a circular velocity of 40 km/s, roughly three times that observed for the dSph’s. Figure 4 shows the “Tully-Fisher” relation for the 11 Galactic satellites compared with the curves predicted from the semi-analytic models employed in Figure 3. The circular velocities are overestimated by a factor of 3–4.

Figure 4. The Tully-Fisher relation for the Galactic satellites (green points) where the velocity dispersions of the spheroidals have been converted to circular velocity assuming isotropic orbits and isothermal potentials. The curves show the Tully-Fisher relation predicted by the semi-analytic models from Figure 3 where we varied the efficiency of feedback to match the numbers of dwarf galaxies.

## III Re-Ionisation

We have just demonstrated that feedback does not solve the overabundance of CDM satellites – clearly some form of stochastic biasing is required. A solution was proposed by Bullock etal (2000) in which only those dark matter halos that have virialised prior to re-ionisation can cool gas and form stars. Once the IGM has been reheated then the smallest CDM halos cannot capture or cool gas and they remain completely dark.

In Figure 5 we mark all the progenitor halos that satisfy the condition for cooling gas prior to z=10, which we will take as the redshift of re-ionisation. We mark particles red if they lie within a region of overdensity larger than 1000. The locations of these particles are subsequently tracked to z=0 and marked in the right panel of Figure 5. Roughly 100 satellites satisfy the density criteria at a redshift z=10 and $`80`$ of these physically merger together to form the very central region of the final galaxy halos. The remaining 20 survive intact and can be found orbiting within the virial radius of the two halos (see Figure 6). The mean radius of the surviving satellites is $`80`$ kpc, which is a factor of 2.5 smaller than the half mass radius of the final halos.

Figure 5. The left panel shows the Local Group simulation at z=10. Marked in red are all those particles that lie in regions with an overdensity larger than 1000. The right panel shows one of the high resolution halos at z=0 and the locations of the red particles marked at z=10.

Figure 6. The smoothed distribution of “starlight” in the Local Group at the present day. I plot only those stars that could form in dark matter halos prior to re-ionisation at z=10. The distribution of these stars is highly biased. Roughly a dozen dark matter dominated satellites orbit within each of the parent halos and they have a spatial distribution that matches the real Local Group. Most of the population II/III stars lie at the very centers of the halos surrounding M31 and the Galaxy. Their half light radius is just a few kiloparsecs (c.f. White & Springel 1999) and their luminosity density falls as $`r^3`$ (c.f. Figure 7).

The final cumulative distribution of satellites within one of the simulated halos is shown in Figure 2 and provides a good match to the corrected observational data points. Several puzzles remain. Why don’t we find any satellites in the Galactic halo with velocity dispersion less than $`7`$ km/s? Is cooling that inefficient below $`10`$ km/s such that we do not find any dark matter dominated systems containing just a handful of stars?

The star formation histories of the Local Group satellites presents a further puzzle. Most of the satellites show evidence for several bursts of star formation, some continuing to the present day. Both re-ionisation and the “essential” feedback have been extremely inefficient at removing gas from these tiny halos that have masses $`10^8M_{}`$.

## IV Resolution issues

Central halo profiles

Have we really converged on the unique central structure of CDM halos? This is a hard numerical calculation to perform since we are always relaxation dominated on small scales. In a hierarchical universe the first halos to collapse will contain just a few particles and have relaxation times much shorter than a Hubble time.

Figure 7 shows the final density profile of one of the high resolution Local Group CDM halos. We also plot the final density profile of those particles that were located within highly non-linear regions at redshifts z=10 and z=20. The central 5 kpc of the halo is dominated by those particles that were in virialised halos at z=10. Most of these halos contained just a few particles and their internal structure is completely dominated by resolution effects. Until we can adequately resolve objects collapsing at z=10, we cannot claim to have converged upon the slope of the density profile at 1–2% of the virial radius.

Figure 7. The solid curve shows the density profile of the high resolution halo shown in Figure 5. The dotted and dashed curves show the density profiles of those particles that lie in regions of overdensity larger than 50 at z=20 and z=10 respectively. The long-dashed curve shows the difference between the solid and short-dashed curves. The radial density profile of the marked particles at z=0 has a gradient of -3 which is similar to that of the Galactic spheroid.

Beam Smearing

Rotation curves of dwarf galaxies first highlighted potential problems with the structure of CDM halos (Moore etal 1994, Flores etal 1994, Burkert 1995). The quality of these data were recently questioned by several authors including van den Bosch & Swaters (2000) used rotation curves from 19 dwarf galaxies to claim that CDM halos are consistent with the data. However, to make this statement these authors had to throw away half of the galaxies and adopt unphysical (zero) mass to light ratios. Furthermore, seven of the remaining nine galaxies require concentration parameters in the range c=3–5 which cannot be obtained in any reasonable $`\mathrm{\Lambda }CDM`$ model. One could rephrase the conclusions of these authors by stating that only 2 galaxies from a sample of 19 are consistent with CDM!

Finally, I show the $`H_\alpha `$ and $`HI`$ rotation curves of the nearby dwarf NGC3109 (Blais-Oullette etal 2001). These data clearly show that beam smearing is not an issue for the nearby dwarf galaxies. Furthermore, only a constant density core can fit these data. CDM profiles with central cusps $`<1`$ are ruled out for any value of the concentration parameter. If CDM is correct then we are forced to conclude that galaxies such as NGC3109, NGC5585, IC2574, etc, are somehow strange and that their disk kinematics are somehow not measuring the mass distribution.

Figure 8. The rotation curve of NGC3109 (Blais-Oullette etal 2000) measured in $`HI`$ (filled circles) and $`H_\alpha `$ (open circles). Beam smearing is clearly not an issue with nearby dwarf galaxies.

Acknowledgments I would like to thank Carlton Baugh for constructing Figure 3 and Figure 4. BM is supported by the Royal Society.

References

Blais-Ouellette, S., Carignan, C. & Amram, P. 2000, AJ, in press, astro-ph/0006449.

Bullock, J.S., Kravtsov, A.V. & Weinberg, D.H. 2000, ApJ, astro-ph/0007295.

Bullock, J.S., Kravstov, A.V. & Weinberg, D.H. 2000, ApJ, 548, 33.

Burkert, A. 1995, Ap.J.Lett., 447, L25.

Cole, S. Lacey, C., Baugh, C. Frenk, C.S. 2000, MNRAS, astro-ph/0007281.

Dubinski, J. & Carlberg, R. 1991, Ap.J., 378, 496.

Flores, R.A. & Primack, J.R. 1994, Ap.J.Lett., 457, L5.

Fukushige, T & Makino, J. 1997, Ap.J.Lett., 477, L9.

Ghigna, S., Moore, B., Governato, F., Lake, G., Quinn, T. & Stadel, J. 1998, M.N.R.A.S., 300, 146.

Jing, Y.P. 2000, ApJ, 535, 30.

Klypin, A., Kravtsov, A.V., Valenzuela, O., & Prada, F. 1999, ApJ, 522, 8.

Moore, B., 1994, Nature, 370, 620.

Moore, B., Ghigna, S., Governato, F., Lake, G., Quinn, T., Stadel, J., & Tozzi, P. 1999, ApJLett, 524, L19.

Moore, B., Quinn, T., Governato, F., Stadel, J., & Lake, G. 1999, MNRAS, 310, 1147.

Moore, B., Calcaneo, C., Quinn, T., Governato, F., Stadel, J., & Lake, G. PhysRevD, submitted.

Navarro, J.F., Frenk, C.S. & White, S.D.M. 1996, Ap.J., 462, 563.

van den Bosch, F.C. & Swaters, R.A. 2000, AJ, submitted.

Warren, S.W., Quinn, P.J., Salmon, J.K. & Zurek, H.W. 1992, Ap.J., 399, 405.

White, S.D.M. 2000, ITP conference, url: online.itp.ucsb.edu/online/galaxy\_c00

White, S.D.M. & Springel, V. 1999, ESO Workshop, eds A. Weiss, T. Abel & V. Hill, astro-ph/9911378.


text:
# Dissipation induced Instabilities and the Mechanical Laser

## Abstract

We study the 1:1 resonance for perturbed Hamiltonian systems with small dissipative and energy injection terms. These perturbations of the 1:1 resonance exhibit dissipation induced instabilities. This mechanism allow us to show that a slightly pumping optical cavity is unstable when one takes into account the dissipative effects. The Maxwell-Bloch equations are the asymptotic normal form that describe this instability when energy is injected through forcing at zero frequency. We display a simple mechanical system, close to the 1:1 resonance, which is a mechanical analog of the Laser.

A fundamental problem in mechanics is the determination of the stability of equilibria. Hamiltonian vector fields can undergo a variety of instabilities as a single bifurcation parameter is varied Abraham . There are two fundamental codimension-one bifurcations of Hamiltonian systems. The first is the stationary or steady state bifurcation, which is characterized by two eigenvalues merging at zero with multiplicity two Marsden . The second is the 1:1 resonance, which is the collision of two pure imaginary eigenvalues (and their complex conjugates) at finite frequencies with multiplicity two LordKelvin . This latter instability has been called by different names depending on the field where it has appeared; for instance, ‘confusion of frequencies’ in mechanical or electrical oscillators Rocard , ‘dispersive instability’ Gibbon  and ‘Hamiltonian Krein-Hopf bifurcation’ Marsden . These instabilities are a consequence of the fact that in the Hamiltonian case, if $`\lambda `$ is an eigenvalue, then LordKelvin  so is $`\lambda `$. The same property is present in time reversible systems, i.e., they exhibit identical generic instabilities. Recently, the instabilities of quasi-reversible systems have also been characterized Clerc1999 ; Clerc1999b , in which the irreversible effects are small and can be considered as perturbative terms close to the instability. In particular, systems close to the quasi-reversible 1:1 resonance are described by the Maxwell-Bloch equations when the energy is injected through a forcing at zero frequency Clerc1999b . The Maxwell-Bloch equations describe the interaction of an electromagnetic field and a collection of two level atoms at an optical cavity Lamb .

The aim of this letter is to study the 1:1 resonance for systems in the neighborhood of a Hamiltonian one; that is, we shall consider perturbed Hamiltonian systems with dissipative and energy injection terms. Near this instability, the dissipative terms are responsible for a spectral bifurcation, i.e., the dissipation induces an instability Marsden1994 . This mechanism allows us to show that a slightly pumping optical cavity is unstable when one takes into account the dissipative effects. The Maxwell-Bloch equations are the asymptotic normal form that describes this instability in the presence of a conserved quantity.We shall display a simple mechanical system, which we call the Mechanical laser, which, close to the 1:1 resonance, is a mechanical analog of the Laser.

For the sake of simplicity, we first consider a Hamiltonian system with two degrees of freedom $`H(p^i,q_i,\left\{\lambda \right\})`$, where $`p^i`$ and $`q_i`$ are canonically conjugate variables, $`\left(i=1,2\right)`$, and $`\left\{\lambda \right\}`$ is a set of parameters. Assume there is a 1:1 resonance for an equilibrium at $`\lambda =\lambda _c`$, i.e., the spectrum at the equilibrium has a pair of pure imaginary eigenvalues of multiplicity two, say at $`\pm i\mathrm{\Omega }`$. Nearby, the instability the system is governed by the normal form Elphick

$`_tA`$ $`=`$ $`i\mathrm{\Omega }A+B`$
$`_tB`$ $`=`$ $`i\mathrm{\Omega }B+f(|A|^2,i\left(AB^{}BA^{}\right),\left\{\lambda \lambda _c\right\})A+`$ (1)
$`ig(|A|^2,i\left(AB^{}BA^{}\right),\left\{\lambda \lambda _c\right\})B,`$
where $`f`$ and $`g`$ are complex polynomial functions. There is a change of variables from the given ones $`\{p^1,q_1,p^2,q_2\}`$ to the new ones (the complex variables $`A,B`$) of the form $`\{p^1,q_1,p^2,q_2\}=A\stackrel{}{\mathrm{\Psi }}+B\stackrel{}{\chi }+\stackrel{}{N}(A,B,A^{},B^{})`$, where $`\stackrel{}{\mathrm{\Psi }}`$ is the eigenvector of the linearized system corresponding to $`i\mathrm{\Omega }`$ and $`\stackrel{}{\chi }`$ is a generalized eigenvector in the Jordan sense, and $`\stackrel{}{N}(A,B,A^{},B^{})`$ is a nonlinear vector. This change of variables is not canonical in general. When one considers rotated variables $`\left(A=e^{i\mathrm{\Omega }t}A^{},B=e^{i\mathrm{\Omega }t}_tA^{}\right)`$ and the dominant terms, the normal form reads (omitting the primes)

$$_{tt}A=\epsilon A+i\delta _tA\alpha |A|^2A$$
(2)
where $`\epsilon `$ is the bifurcation parameter, which is proportional to $`\lambda \lambda _c`$; henceforth we assume $`\epsilon 1`$, $`\delta `$ is the gyroscopic term LordKelvin  or detuning term Lamb  and $`\alpha `$ is an order one parameter. The variables and parameters scale as $`A\epsilon ^{1/2}`$, $`_tA\epsilon `$, $`_{tt}A\epsilon ^{3/2}`$, and $`\delta \epsilon ^{1/2}`$. All terms in this equation are of order $`\epsilon ^{3/2}`$, while the higher order terms are $`\epsilon ^2`$.

When equation (1) is written as a second order system, it has an imaginary term linear in $`A`$, but such a term is not present in (2) because it breaks the eigenvalue symmetry ($`\lambda \lambda `$). The asymptotic normal form (2) has the following Hamiltonian

$$H=_tA_tA^{}\epsilon |A|^2+\frac{\alpha }{2}|A|^4,$$
(3)
with the Poisson-Bracket ($`F`$ and $`G`$ are real)

$$\{F,G\}=\frac{F}{A}\frac{G}{A_t^{}}\frac{G}{A}\frac{F}{A_t^{}}+i\delta \frac{F}{A_t}\frac{G}{A_t^{}}+\mathrm{c}.\mathrm{c}.$$

Thus, to cubic order, the nonlinear change of variables is canonical. The eigenvalues of the zero solution ($`A=0`$) are $`\pm i\left(\delta /2\pm 1/2\sqrt{\delta ^24\epsilon }\right)`$. When $`\epsilon \delta ^2/4`$ is negative, the initial Hamiltonian system has four distinct pure imaginary eigenvalues; $`\epsilon \delta ^2/4`$ equal to zero is the 1:1 resonance at $`\pm \mathrm{\Omega }`$ frequencies and when it is positive, the eigenvalues have nonzero real part. Note that the gyroscopic term is a stabilizing effect LordKelvin .

We now consider this Hamiltonian system under the influence of small dissipative terms. This leads to a new term in the asymptotic normal form as follows:

$$_{tt}A=\epsilon A\left(\mu i\delta \right)_tA\alpha |A|^2A,$$
(4)
where $`\mu `$ is positive and order $`\epsilon ^{1/2}`$.

To study the effects of dissipative terms in (4 ), we consider its characteristic polynomial, which has the form $`\lambda ^4+2\mu \lambda ^3+\left(2\epsilon \mu ^2+\delta ^2\right)\lambda ^2\lambda 2\epsilon \mu +\epsilon ^2`$. Looking for roots of the form $`\lambda =a\pm ib`$ and $`\lambda =c\pm id`$, one recognizes that $`\mu =\left(a+c\right)`$ and $`\epsilon \mu =\left(c\left(a^2+b^2\right)+a\left(c^2+d^2\right)\right)`$. Hence, when $`\epsilon `$ is negative, $`a`$ and $`c`$ are negative, i.e., all eigenvalues are to the left of the imaginary axis. For $`\epsilon `$ positive and $`\epsilon \delta ^2/4`$ negative, the unperturbed system is marginal, but the perturbed one satisfies $`ac<0`$; in addition, the eigenvalues with larger frequency move to the left of the imaginary axis (stable modes) and the others to the right (unstable modes), but the eigenvalues that move furthest away from the axis are stable. Finally, when $`\epsilon >\delta ^2/4`$, the eigenvalues have nonzero real part, and again the stable modes are the furthest from the imaginary axis. Therefore, we are led to consider dissipation induced instabilities Marsden1994 . The destabilizing effects through positive or negative total dissipative perturbation was known a long time ago by Lord Kelvin et al. LordKelvin  From the physical point of view, one can understand this phenomena as follows: when $`\epsilon `$ is negative, the energy (3) has a minimum at the origin, hence when dissipation is added, the solutions near the origin move towards it. By contrast, if $`\epsilon `$ is positive, the energy has a saddle point at the origin and when $`\epsilon <\delta ^2/4`$, this is unstable with an algebraic evolution in time, so when one adds dissipation, this is consistent with the conclusion that the solutions near the origin move away from it, exponentially in time.

Using the preceding analysis, we infer that close to the 1:1 resonance, generically the dissipative terms induce an instability. In the case that the instability happens with null detuning ($`\delta =0`$), the normal form has only a real coefficient and so the dissipative terms do not induce instability. A physical example of this last situation is the weakly dissipative baroclinic instability when the effect of earth’s sphericity is ignored Pedlosky  and another is the Kelvin-Helmholz instability Weissman .

As we shall see in detail later, the Laser is a system that shows dissipation induced instability, but firstly, we need to discuss how the energy is injected to the system, so that the solution that becomes unstable is persistent when the dissipative terms are added. There are two natural ways to inject energy to the modes, namely through forcing at finite frequency or at zero. The latter situation is common in physical systems. In a Hamiltonian system this is only possible if there is a conserved quantity, that is, a zero eigenvalue whose mode is nonlinearly coupled with the other ones. For instance, when there is a cyclic variable, the respective momentum $`Z`$ is conserved. A Hamiltonian system that has a 1:1 resonance in the presence of conserved quantity, leads us to consider an extra equation of the form $`_tZ=0`$ and an additional term $`ZA`$ in the equation ( 4). The system presents different behaviors depending on the value of $`ϵZ`$. When one includes dissipative terms and forcing, the asymptotic normal form reads

$`_{tt}A`$ $`=`$ $`\left(\epsilon Z\right)A\left(\mu i\delta \right)_tA\alpha |A|^2A`$
$`_tZ`$ $`=`$ $`\nu Z+\eta |A|^2`$ (5)
where $`Z\epsilon `$, $`\nu \epsilon ^{1/2}`$ and $`\eta \epsilon ^{1/2}`$. The term $`\eta |A|^2`$ permits a non trivial coupling between the variables. When the unperturbed Hamiltonian system has more modes without resonances between them, the perturbed system is governed by the above equations, since the intensities of the other modes decreases in time. Through a nonlinear change of variables, the previous equations are equivalent to the Maxwell-Bloch equations Clerc1999b . Using a multiscaling method, the dispersive instability with small dissipation is also described by the previous equations Gibbon .

To illustrate how dissipation induced instabilities enter, we consider the semiclassical description of the Laser. This is based on the self-consistent interaction of the electromagnetic field with an active medium within an optical cavity. The electric field is described classically (by the Maxwell equations) and the matter as ensemble of atoms possessing two quantized energy levels; phenomenological terms are added to complete the description. Thus, the system is described by Lamb ; Seigman

$`{\displaystyle \frac{^2E}{t^2}}`$ $`={\displaystyle \frac{^2E}{x^2}}{\displaystyle \frac{^2P}{t^2}}\kappa {\displaystyle \frac{E}{t}},`$
$`{\displaystyle \frac{^2P}{t^2}}`$ $`=\gamma _{}{\displaystyle \frac{P}{t}}(\gamma _{}^2+(1+\delta )^2)PgNE,`$
$`{\displaystyle \frac{N}{t}}`$ $`=\gamma _{}(NN_0)+E\left({\displaystyle \frac{P}{t}}+\gamma _{}P\right),`$ (6)
with periodic boundary condition at the cavity length ($`L`$, $`E(0,t)=E(L,t)`$). Here $`E`$,$`P`$ and $`N`$ are dimensionless quantities, which correspond to linearly polarized electric field, the dipole polarization field and the population inversion. $`\gamma _{}`$, $`\gamma _{}`$ are the decay rate associated to spontaneous emission and interaction between the atoms, $`\kappa `$ is a damping related to the mirror losses, $`\delta `$ is the detuning, $`g`$ is a coupling constant which characterizes the atoms and $`N_0`$ the pump parameter. In the time reversible limit of the above equations, i.e., $`\gamma _{}=\gamma _{}=\kappa =0`$, the system has Hamiltonian density Holm

$$=\left(\frac{D}{4}\frac{1}{\omega }Re\left(a_+a_{}^{}\right)\right)^2+\frac{\left(_xA\right)^2}{2}+|a_+|^2|a_{}|^2,$$

with the Poisson-bracket

$`\{F,G\}`$ $`=`$ $`{\displaystyle \frac{dx}{2}\left(\frac{F}{A}\frac{G}{D}\frac{G}{D}\frac{F}{A}\right)}+`$
$`i{\displaystyle \frac{\omega }{2}}\left({\displaystyle \frac{F}{a_+}}{\displaystyle \frac{G}{a_+^{}}}+{\displaystyle \frac{F}{a_{}}}{\displaystyle \frac{G}{a_{}^{}}}\right)+\mathrm{c}.\mathrm{c}.`$
Where $`N|a_+|^2|a_{}|^2`$, $`P\frac{4}{\omega }Re\left(a_+a_{}^{}\right)`$, $`EA_t`$, and $`\omega =1+\delta `$. The Hamiltonian is just the sum of the electromagnetic energy and the atomic excitation energy Holm .

Changing the cavity length in the time reversible limit of equation (6), leads to a 1:1 resonance for the non-lasing solution ($`E=P=0,`$ $`N=D_0`$), which gives rise to an electromagnetic wave with $`\pm \omega `$ frequencies. Using the slowly varying envelope (WKB) approximation leads to the Maxwell-Bloch equations Lamb .

Figure 1 shows the space-time diagram of the electric field of a numerical simulation of the semiclassical model, plus noise with small intensity, close to the 1:1 resonance without dissipative terms (see fig. 1a) and with dissipative ones (see Fig. 1b). The numerical simulations start with the same initial condition, namely, the no-lasing solution with excited atoms ($`D_0>0`$). It is clear from these pictures that the inclusion of dissipation induces the Laser to respond. One can physically understand what happens, since without dissipation the atoms excited decay through the stimulated emission, i.e., the non-lasing solution becomes unstable very slowly (nonlinear mechanism). Instead, when one takes into account the dissipative terms, the excited atoms decay, for instance through stimulated emission and collisions, i.e., exponentially in time. Note that, one must to pump to $`N_0=D_0`$ in order that the non-lasing solution persists. In brief, with nonzero detuning and slight pumping, the dissipation induces the Laser to respond.

To illustrate the 1:1 resonance in a simple Hamiltonian system, we consider a mechanical system, that we call the Mechanical Laser, which consists of two coupled spherical pendula in a gravitational field, with a support, which can rotate around a vertical axis. The lower pendulum is constrained to move in a plane that is orthogonal to the plane of the upper pendulum (see Fig. 2).

The system rotates with angular velocity $`\dot{\phi }`$ with respect to the vertical. The quantities $`m_1`$, $`m_2`$, $`l_1`$ and $`l_2`$ are the mass and length of the upper and lower pendula respectively, and $`I`$ is the dimensionless moment of inertia of the support. The system will dissipate energy because of friction at the contacts and the motion of the pendulum masses in a fluid (for example, the air) via Stokes’ law. Energy is injected through a constant torque at the upper pendulum pivot point. The governing equations for the angles $`\theta _1\left(t\right)`$ and $`\theta _2\left(t\right)`$ and the vertical angular velocity $`\dot{\phi }`$ read Clerc2001

$$\ddot{\theta }_1=\sigma ^2\mathrm{sin}\theta _1\mathrm{sin}\theta _2\ddot{\theta }_2\sigma ^2\mathrm{sin}\theta _1\mathrm{cos}\theta _2\dot{\theta }_2^2$$

$$2\sigma ^2\mathrm{cos}\theta _1\mathrm{cos}\theta _2\dot{\phi }\dot{\theta }_2+\mathrm{sin}\theta _1\mathrm{cos}\theta _1\dot{\phi }^2$$

$$\sigma ^2\mathrm{cos}\theta _1\mathrm{sin}\theta _2\ddot{\phi }\frac{g}{l}\mathrm{sin}\theta _1\nu _1\dot{\theta }_1,$$

$$\ddot{\theta }_2=\mathrm{sin}\theta _1\mathrm{sin}\theta _2\ddot{\theta }_1\mathrm{cos}\theta _1\mathrm{sin}\theta _2\dot{\theta }_1^2$$

$$+2\mathrm{cos}\theta _1\mathrm{cos}\theta _2\dot{\phi }\dot{\theta }_1+\mathrm{sin}\theta _2\mathrm{cos}\theta _2\dot{\phi }^2$$

$$+\mathrm{sin}\theta _1\mathrm{cos}\theta _2\ddot{\phi }\frac{g}{l}\mathrm{sin}\theta _2\nu _2\dot{\theta }_2,$$

$$\frac{d}{dt}\left\{\begin{array}{c}\left(\mathrm{sin}^2\theta _1+\sigma ^2\mathrm{sin}^2\theta _2\right)\dot{\phi }+\sigma ^2\mathrm{cos}\theta _1\mathrm{sin}\theta _2\dot{\theta }_1\\ \sigma ^2\mathrm{sin}\theta _1\mathrm{cos}\theta _2\dot{\theta }_2+I\dot{\phi }\end{array}\right\}$$

$$=\nu _\phi \left(\dot{\phi }\mathrm{\Omega }\right)\mu _1\mathrm{sin}^2\theta _1\dot{\phi }\mu _2\left(\mathrm{sin}^2\theta _1+\mathrm{sin}^2\theta _2\right)\dot{\phi }.$$

where $`\nu _1`$, $`\nu _2`$, $`\nu _\phi `$, $`\mu _1`$ and $`\mu _2`$ are damping coefficients, $`\sigma =\sqrt{m_2/\left(m_1+m_2\right)}`$ is the relative factor of the energy between the oscillators, and we have written the torque as $`\nu _\phi \mathrm{\Omega }`$. For the sake of simplicity, we have considered the case of pendula of equal lengths ($`l_1=l_2=l`$). When one regards the Hamiltonian limit of the previous equations, the vertical solution or non-lasing solution $`\theta _1=\theta _2=0`$, $`\dot{\phi }=\mathrm{\Omega }_0`$ has a 1:1 resonance when

$$\mathrm{\Omega }_0=\mathrm{\Omega }_c=\sqrt{\frac{g}{l}\frac{\left(m_1+m_2\right)}{m_1}}$$

with frequency $`\omega _c=\pm \sqrt{gm_2/lm_1}`$. The centripetal force is more intense than the gravitational force when $`\mathrm{\Omega }_0\mathrm{\Omega }_g\sqrt{g/l}`$. As a consequence, the Coriolis force exerted by one pendulum on the other, the non-lasing solution is marginal when $`\mathrm{\Omega }_g\mathrm{\Omega }_0\mathrm{\Omega }_c`$. In this region the system is nonlinearly unstable and the system becomes linearly unstable when $`\mathrm{\Omega }_0>\mathrm{\Omega }_c`$; this exhibits a coherence oscillation, which is the signatures of the laser instability.

Near the 1:1 resonance, the coupled pendulum is govern by (5), whereClerc2001

$`\epsilon `$ $`=`$ $`2{\displaystyle \frac{g}{l}}{\displaystyle \frac{\left(\mathrm{\Omega }\mathrm{\Omega }_c\right)}{\mathrm{\Omega }_c}},\text{ }\alpha ={\displaystyle \frac{g}{4l}}\left({\displaystyle \frac{\sigma ^42\sigma ^32\sigma ^2+3}{1\sigma ^2}}\right),`$
$`\delta `$ $`=`$ $`2\sigma \left(\mathrm{\Omega }\mathrm{\Omega }_c\right),\text{ }\nu ={\displaystyle \frac{\nu _\phi }{I}},\text{ }\mu ={\displaystyle \frac{1}{2l^2}}\left(\nu _1+\nu _2\right),`$
$`\eta `$ $`=`$ $`{\displaystyle \frac{1}{I}}\left[\mu _1+\left(1+{\displaystyle \frac{1}{\sigma ^2}}\right)\mu _22{\displaystyle \frac{\nu _\phi }{\sigma ^2I}}\left({\displaystyle \frac{\mathrm{\Omega }}{\mathrm{\Omega }_c}}1\right)\right];`$
where the variables are related to the dominate order by $`A=\left(\sigma \theta _1+i\theta _2\right)\mathrm{exp}(i\omega _ct)/2`$ and $`Z=2\mathrm{\Omega }_c\sigma ^2\left(\mathrm{\Omega }_c\dot{\phi }\right)`$. Transforming equations (5) to the Maxwell-Bloch equations, one obtains the following analogy of the electric field, polarization and population inversion

$`E`$ $`=`$ $`{\displaystyle \frac{\left\{\sigma \theta _1+i\theta _2\right\}}{2}}e^{i\left(\mathrm{\Delta }+\omega _c\right)t}+c.c.`$
$`P`$ $`=`$ $`{\displaystyle \frac{\left\{\sigma \dot{\theta }_1+i\dot{\theta }_2+a\left(\sigma \theta _1+i\theta _2\right)\right\}}{2}}e^{i\left(\mathrm{\Delta }+\omega _c\right)t}+c.c.`$
$`N`$ $`=`$ $`{\displaystyle \frac{2\mathrm{\Omega }_c\sigma ^2}{\alpha }}\left(\dot{\phi }\mathrm{\Omega }_c\right)+\epsilon _0{\displaystyle \frac{\sigma ^2\theta _1^2+\theta _2^2}{4}}`$
where $`\mathrm{\Delta }=\delta \left(\nu \alpha +\eta \right)/\left(2\nu \alpha +2\eta \mu \alpha \right),`$ $`a=\left(\nu \alpha +\eta \right)/\alpha +i2\mathrm{\Delta }`$ and $`\epsilon _0=\left(\epsilon +\mathrm{\Delta }^2\delta \mathrm{\Delta }\right)/\alpha +\left(\nu \alpha +\eta \right)(\mu \alpha \nu \alpha \eta )/\alpha ^3.`$ The mechanical analog of the electric field is a simple function of the pendulum displacements with respect to the vertical. In the Hamiltonian case the polarization is the slow time derivative of the electric field. The mechanical analogy of the stimulated emission mechanism (time reversal effect) is the vertical angular momentum conservation. Thus, when one increases the pendula tilt (the intensity of the electric field enlargement), the angular velocity decreases (the population inversion declines). This is a main ingredient of Laser theory. If the torque or pumping is bigger than the gravitational force ($`\mathrm{\Omega }_0>\mathrm{\Omega }_g`$), i.e., the system is excited in Laser terminology, then the torque that gives rise to an angular velocity slightly higher than $`\mathrm{\Omega }_g`$ leads to a dissipation induced instability.

In summary, we have studied the 1:1 resonance for perturbed Hamiltonian systems with small dissipative and energy injection terms. Nearby, the 1:1 resonance exhibits dissipation induces instabilities. This allow us to show that a slightly pumping optical cavity is unstable when one takes into account the dissipative effects. The Maxwell-Bloch equations are the asymptotic normal form that describe this instability when energy is injected through forcing at zero frequency. We have displayed a simple mechanical system, the Mechanical laser or double spherical pendulum with a support, which, close to the 1:1 resonance, is a mechanical analog of the Laser.

We thank Richard Murray for his encouragement and help in making this work possible and CDS for its hospitality. JEM was partially supported by the National Science Foundation.


text:
# High Energy Neutrino Astronomy: Towards Kilometer-Scale Detectors*footnote **footnote *Talk presented at the International Symposium on High Energy Gamma Ray Astronomy, Heidelberg, June 2000.

## I Introduction

Of all high-energy particles, only neutrinos can directly convey astronomical information from the edge of the universe—and from deep inside the most cataclysmic high-energy processes. Copiously produced in high-energy collisions, travelling at the velocity of light, and not deflected by magnetic fields, neutrinos meet the basic requirements for astronomy. Their unique advantage arises from a fundamental property: they are affected only by the weakest of nature’s forces (but for gravity) and are therefore essentially unabsorbed as they travel cosmological distances between their origin and us.

Many of the outstanding mysteries of astrophysics may be hidden from our sight at all wavelengths of the electromagnetic spectrum because of absorption by matter and radiation between us and the source. For example, the hot dense regions that form the central engines of stars and galaxies are opaque to photons. In other cases, such as supernova remnants, gamma ray bursters, and active galaxies, all of which may involve compact objects or black holes at their cores, the precise origin of the high-energy photons emerging from their surface regions is uncertain. Therefore, data obtained through a variety of observational windows—and especially through direct observations with neutrinos—may be of cardinal importance.

The sun is an intense source of electron neutrinos ($`\nu _e`$), albeit of relatively low energy. Solar neutrino astronomy began with first experiments in the mid-1960s; today there are five complementary neutrino detectors viewing the nuclear reactions in the core of the sun and, at the same time, studying the fundamental properties of neutrinos. The sun remained the only astronomical object studied with neutrinos until neutrinos from Supernova 1987A were observed. These are still the only two sources marking the astronomical neutrino spectrum.

Suggestions to use a large volume of water for high-energy neutrino astronomy were made as early as the 1960s reines . In this case, a muon neutrino ($`\nu _\mu `$) interacts with a hydrogen or oxygen nucleus in the water and produces a muon travelling in nearly the same direction as the neutrino. The blue C̆erenkov light emitted along the muon’s $``$kilometer-long trajectory is detected by strings of photomultiplier tubes deployed deep below the surface. DUMAND, a pioneering project located off the coast of Hawaii, demonstrated that muons could be detected by this technique DUMAND , but the planned detector was never realized. A detector composed of about ninety photomultiplier tubes (or optical modules) located deep in Lake Baikal was the first to demonstrate the detection of neutrino-induced muons in natural water baikal . The European collaborations ANTARES antares  and NESTOR NESTOR  plan to deploy large-area detectors in the Mediterranean Sea within the next few years. The NEMO Collaboration is conducting a site study for a future kilometer-scale detector in the Mediterranean NEMO .

The AMANDA collaboration, situated at the U.S. Amundsen-Scott South Pole Station, has strikingly demonstrated the merits of natural ice as a C̆erenkov detector medium B4 ; albrecht . In 1996, AMANDA was able to observe atmospheric neutrino candidates using only 80 eight-inch photomultiplier tubes B4 . With 1997 data, with 302 optical modules instrumenting approximately 6000 tons of ice, AMANDA extracted several hundred atmospheric neutrino events. AMANDA was the first neutrino telescope with an effective area in excess of 10,000 square meters for TeV muons. In January 2000, AMANDA-II was completed. It consists of 19 string with a total of 677 OMs arranged in concentric circles, with the ten strings from AMANDA forming the central core of the new detector. AMANDA has met the key challenge of neutrino astronomy: it has developed a reliable, expandable, and affordable technology for deploying a kilometer-scale neutrino detector named IceCube.

IceCube and the detectors to be located in the Mediterranean Sea will be complementary in several respects. First and foremost, they will cover different portions of the sky. (From its location at the South Pole, IceCube observes the northern sky.) Because of the different properties of Antarctic ice and ocean water—ice scatters light more but absorbs it less—IceCube should have better energy resolution and a larger effective area for neutrino detection, while the deep sea detectors should have somewhat higher pointing resolution. Finally, one of the key challenges in neutrino astronomy has been the ability to deploy and maintain optical modules far below the surface; here, too, ice and ocean are essentially different.

## II Scientific Goals

The many questions a high-energy neutrino telescope can address naturally begin with astronomy, both nearby in the universe and far away. Fascinating issues in particle physics and other branches of science are accessible as well.

Neutrino telescopes will investigate the engines which power active galaxies, the nature of gamma-ray bursts, and the origin of the highest-energy cosmic rays. They will search for galactic supernovae, for the births of the supermassive black holes which power quasars, and for the annihilation products of halo cold dark matter particles (WIMPS and supersymmetric particles). They can perform coincidence experiments with Earth- and space-based gamma-ray observatories, cosmic ray detectors, and even with future gravitational wave detectors such as LIGO.

For guidance in estimating expected signals, we make use of the observed energy in high energy cosmic-ray protons and nuclei as well as in known sources of non-thermal, high-energy gamma radiation. Some fraction of cosmic rays will interact in their sources, whatever they may be, to produce pions. These interactions may be hadronic collisions with ambient gas or photoproduction with intense photon fields near the high-energy sources. In either case, the neutral pions decay to photons while charged pions include neutrinos among their decay products with spectra related to the observed $`\gamma `$-ray spectra. Estimates based on this relationship show that a kilometer-scale detector is needed to see the neutrino signals snowmass .

The observed fluxes of cosmic rays and gamma rays set the scale of neutrino telescopes. If we, for instance, assume that gamma ray bursts are the cosmic accelerators of the highest energy cosmic rays, one can calculate from textbook particle physics how many neutrinos are produced when the particle beam coexists with the observed MeV photons in the original fireball. We thus predict the observation of order 10–100 neutrinos of PeV energy per year in a detector with a kilometer-square effective area. The rate is somewhat smaller when assuming instead that the observed cosmic ray beam originates near supermassive black holes at the center of active galaxies. Interaction of the cosmic rays with the abundant UV light in the galaxy is the source of the observed neutrinos. For detailed reviews, see reference gaisser .

TeV $`\gamma `$-rays may be produced by radiative processes from accelerated electrons, whereas neutrinos must be produced by hadronic processes. Therefore, high-energy neutrino astronomy has the potential to discriminate between hadronic and electromagnetic models for the intriguing TeV emission from objects such as supernova remnants, gamma-ray burst sources and active galactic nuclei. With this capability, neutrino telescopes may resolve unambiguously one of the oldest puzzles in science: the origin of cosmic rays. TeV photon fluxes determine, independently from the cosmic ray flux, the scale of the neutrino telescope. If one assumes that the observed TeV gamma rays emitted by blazars and supernova remnants are the decay products of neutral pions, one can estimate the accompanying flux of neutrinos in a largely model-independent way. One anticipates neutrino rates of order a few events per year from the Crab and the observed Markarian sources to, possibly, a few hundred from Sgr Aalvarez . The conclusion is, once more, that a kilometer-scale instrument is required to confirm or rule out the hadronic origin of the highest energy photons.

With high-energy neutrino astrophysics, we are poised to open a new window into space and back in time to the highest-energy processes in the universe. The potential scientific payoff of neutrino astronomy arises from the great penetrating power of neutrinos, which allows them to emerge from dense inner regions of energetic sources. It necessarily follows that the expected interaction rates are small, even with a kilometer-scale detector. History has shown, however, that the opening of each new astronomical window has led to unexpected discoveries. Thus, for example, there could be hidden particle accelerators from which only the neutrinos escape. Yet the science these instruments may do which we cannot anticipate, may well fuel the quest for astronomical knowledge in the next century. Unexpected discovery has followed the inauguration of most new astronomical instruments. Large reflecting telescopes on mountaintops led to the discovery of distant galaxies and an expanding universe; radio telescopes found the cosmic microwave background; X-ray and gamma-ray satellites have uncovered a bestiary of awesome cosmological objects. Even the first modest solar-neutrino observatory revealed a paradox about the nature of fundamental interactions which is yet to be fully explained. No one can predict all that a high-energy neutrino observatory will find in the sky—except that it is very likely to amaze us.

## III Large Natural C̆erenkov Detectors

The first generation of neutrino telescopes, launched by the bold decision of the DUMAND collaboration to construct such an instrument, are designed to reach a large telescope area and detection volume for a neutrino threshold of order 10 GeV. This relatively low threshold permits calibration of the novel instrument on the known flux of atmospheric neutrinos. The architecture is optimized for reconstructing the C̆erenkov light front radiated by an up-going, neutrino-induced muon. Only up-going muons made by neutrinos reaching us through the Earth can be successfully detected. The Earth is used as a filter to screen the fatal background of cosmic ray muons. This makes neutrino detection possible over the lower hemisphere of the detector. Up-going muons must be identified in a background of down-going, cosmic ray muons which are more than $`10^5`$ times more frequent for a depth of 1$``$2 kilometers. The method is sketched in Fig. 1.

The optical requirements of the detector medium are severe. A large absorption length is required because it determines the spacings of the optical sensors and, to a significant extent, the cost of the detector. A long scattering length is needed to preserve the geometry of the C̆erenkov pattern. Nature has been kind and offered ice and water as adequate natural C̆erenkov media. Their optical properties are, in fact, complementary. Water and ice have similar attenuation length, with the role of scattering and absorption reversed. Optics seems, at present, to drive the evolution of ice and water detectors in predictable directions: towards very large telescope area in ice exploiting the long absorption length, and towards lower threshold and good muon track reconstruction in water exploiting the long scattering length.

### III.1 Baikal, ANTARES, Nestor and NEMO:Northern Water

Whereas the science is compelling, the real challenge is to develop a reliable, expandable and affordable detector technology. With the termination of the pioneering DUMAND experiment, the efforts in water are, at present, spearheaded by the Baikal experimentbaikal . The Baikal Neutrino Telescope is deployed in Lake Baikal, Siberia, 3.6 km from shore at a depth of 1.1 km. An umbrella-like frame holds 8 strings, each instrumented with 24 pairs of 37-cm diameter QUASAR photomultiplier tubes (PMT). Two PMTs in a pair are switched in coincidence in order to suppress background from natural radioactivity and bioluminescence. Operating with 144 optical modules since April 1997, the NT-200 detector has been completed in April 1998 with 192 optical modules (OM). The Baikal detector is well understood, and the first atmospheric neutrinos have been identified.

The Baikal site is competitive with deep oceans, although the smaller absorption length of C̆erenkov light in lake water requires a somewhat denser spacing of the OMs. This does, however, result in a lower threshold which may be a definite advantage, for instance for oscillation measurements and WIMP searches. They have shown that their shallow depth of 1 kilometer does not represent a serious drawback. By far the most significant advantage is the site with a seasonal ice cover which allows reliable and inexpensive deployment and repair of detector elements from a stable platform.

With data taken with 96 OMs only, they have shown that atmospheric muons can be reconstructed with sufficient accuracy to identify atmospheric neutrinos; see Fig. 2. The neutrino events are isolated from the cosmic ray muon background by imposing a restriction on the chi-square of the C̆erenkov fit, and by requiring consistency between the reconstructed trajectory and the spatial locations of the OMs reporting signals. In order to guarantee a minimum lever arm for track fitting, they only consider events with a projection of the most distant channels on the track larger than 35 meters. This does, of course, result in a higher threshold.

In the following years, NT-200 will be operated as a neutrino telescope with an effective area between $`10^35\times 10^3`$ m2, depending on energy. Presumably too small to detect neutrinos from extraterrestrial sources, NT-200 will serve as the prototype for a larger telescope. For instance, with 2000 OMs, a threshold of $`1020`$ GeV and an effective area of $`5\times 10^410^5`$ m2, an expanded Baikal telescope would fill the gap between present underground detectors and planned high threshold detectors of cubic kilometer size. Its key advantage would be low threshold.

The Baikal experiment represents a proof of concept for deep ocean projects. These have the advantage of larger depth and optically superior water. Their challenge is to find reliable and affordable solutions to a variety of technological challenges for deploying a deep underwater detector. Several groups are confronting the problem; both NESTOR and ANTARES are developing rather different detector concepts in the Mediterranean.

The NESTOR collaborationNESTOR , as part of a series of ongoing technology tests, is testing the umbrella structure which will hold the OMs. They have already deployed two aluminum “floors”, 34 m in diameter, to a depth of 2600 m. Mechanical robustness was demonstrated by towing the structure, submerged below 2000 m, from shore to the site and back. These tests should soon be repeated with fully instrumented floors. The actual detector will consist of a tower of 12 six-legged floors vertically separated by 30 m. Each floor contains 14 OMs with four times the photocathode area of the commercial 8 inch photomultipliers used by AMANDA and ANTARES.

The detector concept is patterned along the Baikal design. The symmetric up/down orientation of the OMs will result in uniform angular acceptance and the relatively close spacings in a low threshold. NESTOR does have the advantage of a superb site off the coast of Southern Greece, possibly the best in the Mediterranean. The detector can be deployed below 3.5 km relatively close to shore. With the attenuation length peaking at 55 m near 470 nm the site is optically superior to that of all other deep water sites investigated for neutrino astronomy.

The ANTARES collaborationantares  is investigating the suitability of a 2400 m-deep Mediterranean site off Toulon, France. The site is a trade-off between acceptable optical properties of the water and easy access to ocean technology. Their detector concept indeed requires remotely operated vehicles for making underwater connections. First results on water quality are very encouraging with an attenuation length of 40 m at 467 nm and a scattering length exceeding 100 m. Random noise exceeding 50 khz per OM is eliminated by requiring coincidences between neighboring OMs, as is done in the Lake Baikal design. Unlike other water experiments, they will point all photomultipliers sideways in order to avoid the effects of biofouling. The problem is significant at the Toulon site, but only affects the upper pole region of the OM. Relatively weak intensity and long duration bioluminescence results in an acceptable deadtime of the detector. They have demonstrated their capability to deploy and retrieve a string, and have reconstructed down-going muons with 8 OMs deployed on the test string.

With the study of atmospheric neutrino oscillations as a top priority, they had planned to deploy in 2001-2003 10 strings instrumented over 400 m with 100 OMs. After study of the underwater currents they decided that they can space the strings by 100 m, and possibly by 60 m. The ANTARES detector will consist of 13 strings, each equipped with 30 storeys and 3 PMTs per storey. The large photocathode density of the array will allow the study of atmospheric neutrino oscillations in the range $`255<L/E<2550\mathrm{km}\mathrm{GeV}^1`$ with neutrinos in the energy range $`5<E_\nu <50`$ GeV. This detector will have an area of about $`3\times 10^4\mathrm{m}^2`$ for 1 TeV muons — similar to AMANDA-II — and is planned to be fully deployed by the end of 2003.

A new R&D initiative based in Catania, Sicily has been mapping Mediterranean sites, studying mechanical structures and low power electronics. One must hope that with a successful pioneering neutrino detector of $`10^3\mathrm{km}^3`$ in Lake Baikal, a forthcoming $`10^2\mathrm{km}^3`$ detector near Toulon, the Mediterranean efforts will converge on a $`10^1\mathrm{km}^3`$ detector possibly at the NESTOR sitespiro . For neutrino astronomy to become a viable science several of these, or other, projects will have to succeed besides AMANDA. Astronomy, whether in the optical or in any other wave-band, thrives on a diversity of complementary instruments, not on “a single best instrument”. When, for instance, the Soviet government tried out the latter method by creating a national large mirror project, it virtually annihilated the field.

### III.2 AMANDA: Southern Ice

Construction of the first-generation AMANDA detector was completed in the austral summer 96–97. It consists of 302 optical modules deployed at a depth of 1500–2000 m; see Fig. 3. Here the optical module consists of an 8-inch photomultiplier tube and nothing else. It is connected to the surface by a cable which transmits the HV as well as the anode current of a triggered photomultiplier. The instrumented volume and the effective telescope area of this instrument matches those of the ultimate DUMAND Octagon detector which, unfortunately, could not be completed.

Depending on depth, the absorption length of blue and UV light in the ice varies between 85 and 225 metres. The effective scattering length, which combines the mean-free path $`\lambda `$ with the average scattering angle $`\theta `$ as $`\frac{\lambda }{(1cos\theta )}`$, varies from 15 to 40 metres science . Because the absorption length of light in the ice is very long and the scattering length relatively short, many photons are delayed by scattering. In order to reconstruct the muon track one uses maximum likelihood methods, which take into account the scattering and absorption of photons as determined from calibration measurements B4 . A Bayesian formulation of the likelihood, which takes into account the much larger rate of down-going tracks relative to up-going signal, has been particularly effective in decreasing the chance for a down-going muon to be misreconstructed as upward-going.

Other types of events that might appear to be up-going muons must also be considered and eliminated. Rare cases, such as muons which undergo catastrophic energy loss, for instance through bremsstrahlung, or that are coincident with other muons, must be investigated. To this end, a series of requirements or quality criteria, based on the characteristic time and spatial pattern of photons associated with a muon track and the response of the detector, are applied to all events that, in the first analysis, appear to be up-going muons. For example, an event which has a large number of optical modules hit by photons unscattered relative to the expected Cherenkov times of the reconstructed track, has a high quality. By making these requirements (or “cuts”) increasingly selective we eliminate correspondingly more of the background of false up-going events while still retaining a significant fraction of the true up-going muons, i.e., the neutrino signal. Two different and independent analyses of the same set of data covering 138 days of observation have been undertaken. These analyses yielded comparable numbers of upgoing muons (153 in analysis A, 188 in analysis B). Comparison of these results with their respective Monte Carlo simulations shows that they are consistent with each other in terms of the numbers of events, the number of events in common, and the expected properties of atmospheric neutrinos.

In Fig. 4, from analysis A, the experimental events are compared to simulations of background and signal as a function of the (identical) quality requirements placed on the three types of events: experimental data, simulated up-going muons from atmospheric neutrinos, and a simulated background of down-going cosmic-ray muons. For simplicity in presentation, the levels of the individual types of cuts have been combined into a single parameter representing the overall event quality, and the comparison is made in the form of ratios. Figure 4 shows events for which the quality level is 4 and higher. As the quality level is increased further, the ratios of simulated background to experimental data and experimental data to simulated signal both continue their rapid decrease, the former toward zero and the latter toward unity. Over the same range, the ratio of experimental data to the simulated sum of background and signal remains near unity. At an event quality of 6.9 there are 153 events in the sample of experimental data and the ratio to predicted signal is 0.7. The conclusion is that the quality requirements have reduced the events from misreconstructed down-going muons in the experimental data to a negligible fraction of the signal and that the experimental data behave in the same way as the simulated atmospheric neutrino signal for events that pass the stringent cuts. The remaining instrumental background is estimated at $`15\pm 7`$ percent of the signal.

The estimated uncertainty on the number of events predicted by the signal Monte Carlo simulation, which includes uncertainties in the high-energy atmospheric neutrino flux, the in-situ sensitivity of the optical modules, and the precise optical properties of the ice, is +40% $``$50%. The observed ratio of experiment to simulation (0.7) and the expectation (1.0) therefore agree within errors.

The shape of zenith angle distribution from analysis B is compared to a simulation of the atmospheric neutrino signal in Fig. 5 in which the two distributions have been normalized to each other. The variation of the measured rate with zenith angle is reproduced by the simulation to within the statistical uncertainty. Note that the tall geometry of the detector strongly influences the dependence on zenith angle in favor of more vertical muons.

Estimates of the energies of the up-going muons (based on simulations of the number of optical modules that participate in an event) indicate that the energies of these muons are in the range from 100 GeV to $`1\mathrm{TeV}`$. This is consistent with their atmospheric neutrino origin.

The agreement between simulation and experiment shown in Figures 4 and 5 taken together with other comparisons of measured and simulated eventsnature00  leads to the conclusion that the up-going muon events observed by AMANDA are produced mainly by atmospheric neutrinos.

The arrival directions of the neutrinos observed in both analyses are shown in Fig. 6. A statistical analysis indicates no evidence for point sources in this sample. An estimate of the energies of the up-going muons (based on simulations of the number of reporting optical modules) indicates that all events have energies consistent with an atmospheric neutrino origin. This enables us to reach a level of sensitivity to a diffuse flux of high energy extra-terrestrial neutrinos of order $`dN/dE_\nu =10^6E_\nu ^2\mathrm{cm}^2\mathrm{s}^1\mathrm{sr}^1\mathrm{GeV}^1,`$ assuming an $`E^2`$ spectrum. At this level we would exclude a variety of theoretical models which assume the hadronic origin of TeV photons from active galaxies and blazars. Searches for neutrinos from gamma-ray bursts, for magnetic monopoles, and for a cold dark matter signal from the center of the Earth are also in progress and, with only 138 days of data, yield limits comparable to or better than those from smaller underground neutrino detectors that have operated for a much longer period.

Data are being taken now with the larger array, AMANDA-II. New surface electronics consolidates several triggering functions and adds functionality. New scalers were installed that provide millisecond resolution — important for Supernova studies. Several technologies were deployed to evaluate their utility and readiness for future expansion to larger systems.

Yet, the fluxes of very high energy neutrinos predicted by theoretical models or derived from the observed flux of ultra high energy cosmic rays are sufficiently low that a neutrino detector having an effective area up to a square kilometer is required for their observation and studysnowmass . Plans are therefore being made for a much larger detector, IceCube, consisting of 4800 photomultipliers to be deployed on 80 strings. This proposed neutrino telescope would have an effective area of $`1\mathrm{km}^2`$, an energy threshold near 100 GeV and pointing accuracy for muons better than one degree for high energy events. In conclusion, the observation of atmospheric neutrinos reported here for AMANDA is a significant step toward establishing the field of neutrino astronomy first envisioned over 40 years ago.

Although neutrino telescopes are intended primarily for TeV (and higher) energy neutrino astronomy, AMANDA and IceCube also have the potential to detect the burst of several-MeV neutrinos from a galactic supernova and possibly even an extragalactic burst of low-energy neutrinos associated with the birth of a massive black hole. By continuously monitoring the summed counting rate of its 4800 optical modules, ice detectors will monitor the sky for such cataclysmic phenomena. The interactions of this host of low-energy neutrinos will be distributed uniformly throughout the detector, a signal rising above the background noise level. The very low noise rate for an optical module in ice (as opposed to the very high noise rates in ocean water) makes it possible to continue this “supernova watch,” already begun by AMANDA ralf .

## Acknowledgments

This review has been made possible by contributions from Jaime Alvarez, Serap Tilav, and my AMANDA collaborators. This research was supported in part by the U.S. Department of Energy under Grant No. DE-FG02-95ER40896 and in part by the University of Wisconsin Research Committee with funds granted by the Wisconsin Alumni Research Foundation.