text:
# Untitled Document

Copenhagen Computation:

How I Learned to Stop Worrying and Love Bohr

N. David Mermin

Laboratory of Atomic and Solid State Physics

Cornell University, Ithaca, NY 14853-2501

To celebrate the 60th birthday of Charles H. Bennett I (1) publicly announce my referee reports for the original dense coding and teleportation papers, (2) present a very economical solution to the Bernstein-Vazirani problem that does not even hint at interference between multiple universes, and (3) describe how I inadvertently reinvented the Copenhagen interpretation in the course of constructing a simple, straightforward, and transparent introduction to quantum mechanics for computer scientists.

1. Preface: present at the birth

David DiVincenzo, Patrick Hayden, and Barbara Terhal  have designated me the “midwife of teleportation” in recognition of my having written a favorable referee’s report on the discovery paper  and having advised the editors that the proposed terminology made sense. Though this honorific raises vexing biological questions — can something with six fathers and no mother be brought forth by a midwife? — I accept the title with pride. As midwife it seemed appropriate for me to read my referee’s report at the Bennett 60th Birthday Symposium, attended, as it was, by all six fathers. I reproduce it here too, since it shows me to have had a taste for Copenhagen Computation (about which more below) even before Chris Fuchs  got to work on me.

Referee’s Report: Bennett et al., "Teleporting$`\mathrm{}`$" LZ4539
This is a charming, readable, thought-provoking paper.
It presents a novel application of EPR correlations. The
character of the quantum state (how much is inherent in
the physical system, how much is a representation of our
knowledge) is still an extremely elusive notion. This
novel method for duplicating a quantum state somewhere
else by a combination of quantum correlations and classical
information will become an important one of the
intellectual tools available to anybody trying to clear
up this murkiness.
While hunting down the above report I discovered, to my amazement that the year before I had also refereed the discovery paper on dense coding . (I was under the impression that I had paid no attention whatever to dense coding until shortly before its deconstruction in 2002 .)

Bennett and Wiesner, "Communication via one-and
two-particle$`\mathrm{}`$" LT4749
Your question was: Does this qualify as "strikingly
different" enough to publish? I have never read anything
like it, and I have read a lot on EPR, though far from
everything ever written. So as far as I know it is
different.
But strikingly? The argument is very simple, so
shouldn’t the point be obvious? After reading the paper
I put it aside and spent the next week working hard on
something totally unrelated. Every now and then I would
introspect to see if some way of looking at the argument
had germinated that reduced it to a triviality. None had.
Last night I woke up at 3am, fascinated and obsessed with
it. Couldn’t get back to sleep. That’s my definition of
"striking".
So I say it’s strikingly different and I say publish it.
Rereading these old reports reminded me that the myth that referees relish their power to reject papers is off the mark. Writing a favorable report for a good paper is sheer pleasure. Negative reports are no fun at all.

2. Prologue: Bernstein-Vazirani without parallel universes

The Bernstein-Vazirani problem presents one with a black-boxed subroutine, shown in Figure 1, whose action on $`n+1`$ qubits is that of a unitary transformation $`\text{U}_a`$ which takes the computational basis state $`|x_n|y_1`$ of an $`n`$-Qbit  input register and 1-Qbit output register into the state $`|x_n|yxa_1`$. Here $``$ denotes addition modulo 2, $`xa`$ denotes the bitwise modulo 2 inner product of the two $`n`$-bit numbers $`x`$ and $`a`$ ($`xa=x_{n1}a_{n1}\mathrm{}x_1a_1x_0a_0`$), and $`a`$ is some fixed but unknown $`n`$-bit integer with binary expansion $`a=a_{n1}\mathrm{}a_1a_0`$. The problem is to find the smallest number of invocations of the black box needed to learn $`a`$.

Figure 1. The black-boxed Bernstein-Vazirani subroutine $`\text{U}_a`$. The heavy wire is the $`n`$-Qbit input register; the light wire is the 1-Qbit output register.

If the subroutine is applied to $`x=2^j`$, the output register will be flipped if and only if $`a_j=1`$, so a classical computer can determine $`a`$ with $`n`$ calls of the subroutine. Evidently there is no classical way to learn $`a`$ with fewer than $`n`$ calls, since one needs $`n`$ independent linear relations among the bits of $`a`$. But with a quantum computer one can find $`a`$ with just a single call of the subroutine, whatever the size of $`n`$.

This remarkable trick is done by applying a Hadamard transformation,

$$\text{H}|\mathrm{\hspace{0.17em}0}=\frac{1}{\sqrt{2}}\left(|\mathrm{\hspace{0.17em}0}+|\mathrm{\hspace{0.17em}1}\right),\text{H}|\mathrm{\hspace{0.17em}1}=\frac{1}{\sqrt{2}}\left(|\mathrm{\hspace{0.17em}0}|\mathrm{\hspace{0.17em}1}\right),$$
$`(1)`$
to every one of the $`n+1`$ Qbits both before and after the application of $`\text{U}_a`$, as shown in Figure 2. If one initializes the input register to the state $`|\mathrm{\hspace{0.17em}0}_n`$ and the output register to the state $`|\mathrm{\hspace{0.17em}1}`$, then at the end of this process the input register is guaranteed to be in the $`n`$-Qubit state $`|a_n`$. So $`a`$ can be learned by measuring each Qbit of the output register in the computational basis.

Figure 2. Quantum solution of the Bernstein-Vazirani problem. $`\text{H}_n`$ is an $`n`$-fold tensor product of 1-Qbit Hadamards.

The conventional explanation for why this works goes like this:

(i) Applying Hadamards to every Qbit of an input register initially in the $`n`$-Qbit state $`|\mathrm{\hspace{0.17em}0}_n`$ results in a uniform superposition of all possible inputs:

$$\text{H}_n|\mathrm{\hspace{0.17em}0}_n=2^{n/2}\underset{0x<2^n}{}|x_n.$$
$`(2)`$
(ii) Preparing the output register in the state $`\text{H}|\mathrm{\hspace{0.17em}1}`$ converts a bit-flip into a change of phase (specifically a multiplication by $`1`$).

(iii) Another application of Hadamards to the input register after the application of $`\text{U}_a`$ introduces additional $`x`$-dependent phases according to the rule

$$\text{H}_n|x_n=2^{n/2}\underset{0z<2^n}{}(1)^{xz}|z_n.$$
$`(3)`$
(iv) A little arithmetic now reveals that the combined phases lead to complete destructive interference for every term characterizing the input register in the final superposition except for the single state $`|a_n`$.

This process is usually described as an application of massive quantum parallelism followed by destructive interference among all the unfavorable outcomes. People with overactive imaginations like to say that step (i) initializes a computer in each of $`2^n`$ parallel universes to each of the $`2^n`$ possible inputs. The remaining steps are cunningly designed to produce destructive interference among all those $`2^n`$ universes, in just such a way as to lead in every single universe to the presence of $`a`$ in the input register at the end of the process.

There is, however, a much simpler way to understand why the circuit in Figure 2 behaves as advertised, which offers no hint of this metaphysical extravaganza. This approach merely notes that the effect of Hadamards on the basic 2-Qbit controlled-NOT (cNOT) gate, defined in Figure 3, is just to interchange the control and target Qbits, as shown in Figure 4. This follows from the fact that $`\text{H}^2=\text{1}`$ and $`\text{H}\text{X}\text{H}=\text{Z}`$, where

$$\text{X}|\mathrm{\hspace{0.17em}0}=|\mathrm{\hspace{0.17em}1},\text{X}|\mathrm{\hspace{0.17em}1}=|\mathrm{\hspace{0.17em}0},\text{Z}|\mathrm{\hspace{0.17em}0}=|\mathrm{\hspace{0.17em}0},\text{Z}|\mathrm{\hspace{0.17em}1}=|\mathrm{\hspace{0.17em}1},$$
$`(4)`$
and the fact that controlled-Z is symmetric under interchange of target and control Qbits.

Figure 3. The 2-Qbit cNOT gate.

Figure 4. Hadamards interchange target and control Qbits.

The action of $`\text{U}_a`$ shown in Figure 1 is identical to the action of a collection of cNOT gates — one for each non-zero bit of $`a`$. They all target the output register, and are controlled by just those Qbits representing bits of $`x`$ that correspond to non-zero bits of $`a`$. This is illustrated in Figure 5 for $`n=5`$ and $`a=11010`$. Since sandwiching cNOT gates between Hadamards interchanges the control and target Qbits and since H is its own inverse, the magic of Bernstein-Vazirani follows at once, as shown in Figure 6.

Figure 5. The black-boxed Bernstein-Vazirani oracle (shown for the case $`n=5`$, $`a=11010`$) behaves as if it contained a collection of cNOT gates.

Figure 6. Because $`\text{U}_a`$ behaves as a collection of cNOT gates, because Hadamards reverse the action of cNOT gates, and because the output register has been set to 1, we have a simple explanation, requiring only one universe, for why $`a`$ can be determined with only one invocation of $`\text{U}_a`$.

Notice that the only way quantum mechanics enters is through the ability of Hadamards to reverse the action of cNOT gates, as illustrated in Figure 4. Since the same trick can be done classically with 2-Qbit SWAP gates, as shown in Figure 7, the magic of quantum mechanics here lies entirely in the possibility it offers for reversing the roles of target and control Qbits using only 1-Qbit local operations. (If the six pairs of vertically separated Hadamards in the middle circuit of Figure 6 were vertically linked into irreducibly 2-Qbit gates, then they could no longer be moved to the extreme right and left of the circuit without leaving any traces in the central part, as in the upper half of Figure 6.)

Figure 7. Classical SWAP gates also invert the action of a cNOT gate just as Hadamards do (Figure 4). The power of Hadamards over SWAPs is that the Hadamards can do it without requiring any interaction between the two Qbits. If swap gates are used to invert the cNOT gates in the final form of Figure 6, they get all tangled up with each other when one attempts move them to the edges of the figure, as one can do with the Hadamards.

3. How I invented the Copenhagen interpretation while teaching quantum mechanics to computer scientists.

\[I\]n our description of nature the purpose
is not to disclose the real essence of the phenomena
but only to track down, so far as it is possible,
relations between the manifold aspects of our experience.

— Niels Bohr

For the past few years I have taught a course in quantum computation suitable for computer scientists having no background in physics . My first challenge was to develop a minimalist introduction to quantum mechanics that straightforwardly conveyed in a few lectures everything a mathematically sophisticated student needed to know to understand discussions like, for example, that of the preceding section.

The advantage of teaching an approach to a subject as you develop it is that you get striking demonstrations of the ways in which it does and doesn’t work. After several iterations and countless revisions, reorganizations, and reconstructions, the process seemed to be converging. It was only then that I realized that the unproblematic, no-nonsense, lucid, practical pedagogical approach that had so painfully evolved out of my clumsy initial attempts, was nothing but the standard Copenhagen interpretation. What follows, therefore, is my vision of why quantum computation, far from demonstrating the existence of the multiverse, provides the simplest and most compelling example of a major application of quantum mechanics which the Copenhagen point of view fits like a glove.

We begin with a silly formulation of ordinary non-quantum classical computation, based on representing the integers less than $`N`$ as orthonormal vectors in $`N`$ dimensions:

$$0\left(\begin{array}{c}1\\ 0\\ 0\\ 0\\ 0\\ \mathrm{}\end{array}\right),1\left(\begin{array}{c}0\\ 1\\ 0\\ 0\\ 0\\ \mathrm{}\end{array}\right),2\left(\begin{array}{c}0\\ 0\\ 1\\ 0\\ 0\\ \mathrm{}\end{array}\right),3\left(\begin{array}{c}0\\ 0\\ 0\\ 1\\ 0\\ \mathrm{}\end{array}\right),4\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ 1\\ \mathrm{}\end{array}\right),\mathrm{}.$$
$`(5)`$
This clumsy form takes on a rather simpler structure if $`N`$ is a power of 2, so we specialize to the case $`N=2^n`$. When $`n=1`$ we have only two such vectors, which we denote by a more compact pair of symbols due to Dirac:

$$\left(\begin{array}{c}1\\ 0\end{array}\right)=|\mathrm{\hspace{0.17em}0},\left(\begin{array}{c}0\\ 1\end{array}\right)=|\mathrm{\hspace{0.17em}1}.$$
$`(6)`$
To manipulate these two numbers in a computer it is necessary to represent then by a physical system having two distinguishable configurations. Continuing to follow Dirac, we call any such physical system a Cbit  (“C” for “classical”). The vectors $`|\mathrm{\hspace{0.17em}0}`$ and $`|\mathrm{\hspace{0.17em}1}`$ associated with these two configurations are called the states of the Cbit.

If we have two Cbits ($`n=2`$) their four states conveniently decompose into the tensor product of two 1-Cbit states:

$$|\mathrm{\hspace{0.17em}0}_2=\left(\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right)=\left(\begin{array}{c}1\\ 0\end{array}\right)\left(\begin{array}{c}1\\ 0\end{array}\right)=|\mathrm{\hspace{0.17em}0}|\mathrm{\hspace{0.17em}0}=|\mathrm{\hspace{0.17em}00},$$

$$|\mathrm{\hspace{0.17em}1}_2=\left(\begin{array}{c}0\\ 1\\ 0\\ 0\end{array}\right)=\left(\begin{array}{c}1\\ 0\end{array}\right)\left(\begin{array}{c}0\\ 1\end{array}\right)=|\mathrm{\hspace{0.17em}0}|\mathrm{\hspace{0.17em}1}=|\mathrm{\hspace{0.17em}01},$$

$$|\mathrm{\hspace{0.17em}2}_2=\left(\begin{array}{c}0\\ 0\\ 1\\ 0\end{array}\right)=\left(\begin{array}{c}0\\ 1\end{array}\right)\left(\begin{array}{c}1\\ 0\end{array}\right)=|\mathrm{\hspace{0.17em}1}|\mathrm{\hspace{0.17em}0}=|\mathrm{\hspace{0.17em}10},$$

$$|\mathrm{\hspace{0.17em}3}_2=\left(\begin{array}{c}0\\ 0\\ 0\\ 1\end{array}\right)=\left(\begin{array}{c}0\\ 1\end{array}\right)\left(\begin{array}{c}0\\ 1\end{array}\right)=|\mathrm{\hspace{0.17em}1}|\mathrm{\hspace{0.17em}1}=|\mathrm{\hspace{0.17em}11}.$$
$`(7)`$
The last two forms in each line provide some simpler notations for these 2-Cbit states. Pause to admire how the quantum mechanical practice of representing the states of composite systems by the tensor product of the subsystem states emerges automatically from the trivial representation of integers introduced in (5).

The tensor product extends straightforwardly to many Cbits: states of $`n`$ Cbits can be expressed as tensor products of $`n`$ 1-Cbit states. For example,

$$|\mathrm{\hspace{0.17em}5}_3=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ 1\\ 0\\ 0\end{array}\right)=\left(\begin{array}{c}0\\ 1\end{array}\right)\left(\begin{array}{c}1\\ 0\end{array}\right)\left(\begin{array}{c}0\\ 1\end{array}\right)=|\mathrm{\hspace{0.17em}1}|\mathrm{\hspace{0.17em}0}|\mathrm{\hspace{0.17em}1}=|\mathrm{\hspace{0.17em}101}.$$
$`(8)`$
While the operation X defined in (4) makes perfect sense for Cbits (representing the logical NOT) the operation Z makes no sense at all, since we have assigned no meaning to the sign of the state-vector that describes a Cbit. Nevertheless, combinations of operators Z on pairs of Cbits can be classically meaningful. For example

$$\frac{1}{2}\left(\text{1}+\text{Z}\text{Z}\right)\mathrm{projects}\mathrm{on}|\mathrm{\hspace{0.17em}0}|\mathrm{\hspace{0.17em}0},|\mathrm{\hspace{0.17em}1}|\mathrm{\hspace{0.17em}1},$$

$$\frac{1}{2}\left(\text{1}\text{Z}\text{Z}\right)\mathrm{projects}\mathrm{on}|\mathrm{\hspace{0.17em}0}|\mathrm{\hspace{0.17em}1},|\mathrm{\hspace{0.17em}1}|\mathrm{\hspace{0.17em}0}.$$
$`(9)`$
This leads directly to a representation of the SWAP operator S that takes the 2-Cbit state $`|x|y`$ to $`|y|x`$:

$$\begin{array}{cc}\hfill \text{S}=& \frac{1}{2}\left(\text{1}+\text{Z}\text{Z}\right)+\left(\text{X}\text{X}\right)\frac{1}{2}\left(\text{1}\text{Z}\text{Z}\right)\hfill \\ \hfill =& \frac{1}{2}\left(\text{1}+\text{Z}\text{Z}+\text{X}\text{X}\text{Y}\text{Y}\right),\text{Y}=\text{X}\text{Z}.\hfill \end{array}$$
$`(10)`$
Pause also to admire the simplicity of this classical derivation of the form of the quantum-mechanical exchange operator, compared with the standard derivation based on angular momentum technology. Note also the further simplicity introduced into (10) by incorporating an additional factor of $`i`$ into the definition of Y (which also makes it hermitian, like X and Z.) With examples like this one can motivate the utility of extending the notion of states to include multiplication by complex scalars, leading to the generalization from Cbits to Qbits.

Qbits are physical systems characterized by states which fully exploit the entire $`2^n`$ dimensional complex vector space spanned by the $`2^n`$ orthonormal Cbit states. Nature has been kind enough to present us with many examples of them. The general state $`|\mathrm{\Psi }`$ of $`n`$ Qbits is any unit vector:

$$|\mathrm{\Psi }=\underset{0x<2^n}{}a_x|x_n,\underset{0x<2^n}{}|a_x|^2=1.$$
$`(11)`$
With one (extremely important) exception all operations on Qbits are reversible. Since the exception (“measurement”) has no nontrivial analog for Cbits, in comparing Qbits and Cbits it suffices to consider only reversible operations on Cbits. The only reversible operations on the $`2^n`$ Cbit states are their $`(2^n)!`$ possible permutations. But the general operation nature allows us to perform on $`n`$-Qbit states is any linear norm-preserving transformation,

$$|\mathrm{\Psi }\text{U}|\mathrm{\Psi },\text{U}\mathrm{unitary},$$
$`(12)`$
as shown schematically in Figure 8.

Figure 8. Circuit diagram representing a unitary gate U acting on $`n`$ Qbits.

While Qbits are far more versatile than Cbits in their range of states and the operations one can perform on them, the usefulness of their versatility is highly constrained by one important difference between Qbits and Cbits. Learning the state $`|x_n`$ of $`n`$ Cbits is unproblematic: you just look to see which of the $`2^n`$ possible states $`|x_n`$ it is. In contrast, learning the state $`|\mathrm{\Psi }_n=_xa_x|x_n`$ associated with $`n`$ Qbits is impossible. Given the Qbits there is nothing you can do to them to reveal their state.

To extract any information from Qbits one must “make a measurement”. This consists of sending the Qbits through a “measurement gate”. If the state of the Qbits is (11) then the measurement gate signals $`x`$ with probability $`p=|a_x|^2`$. After $`x`$ is signaled, the state associated with the Qbits must be taken to be $`|x_n`$. These features of an $`n`$-Qbit measurement gate and the fact that it represents the only way to extract information from the $`n`$-Qbits is known as the Born rule. The Born rule is illustrated in Figure 9.

Figure 9. Circuit diagram representing an $`n`$-Qbit measurement gate.

As defined above “measurement” is always in the computational basis. One loses nothing by this simplification, since measurement in any other basis can be described as measurement in the computational basis, preceded by an appropriate unitary transformation. What one gains, pedagogically, is the need to invoke only a single variety of measurement gate and, as noted below, only 1-Qbit measurement gates.

A somewhat stronger version of the Born rule plays a crucial role in quantum computation, though it is rarely explicitly mentioned in most standard quantum mechanics texts. The stronger form applies when one measures only a single one of $`n`$ Qbits. The state $`|\mathrm{\Psi }`$ of all $`n`$ Qbits can always be represented in the form

$$|\mathrm{\Psi }=a_0|\mathrm{\hspace{0.17em}0}|\mathrm{\Phi }_0+a_1|\mathrm{\hspace{0.17em}1}|\mathrm{\Phi }_1,|a_0|^2+|a_1|^2=1,$$
$`(13)`$
where the Qbit to be measured appears on the left and where $`|\mathrm{\Phi }_0`$ and $`|\mathrm{\Phi }_1`$ are normalized but not necessarily orthogonal states of the $`n1`$ unmeasured Qbits. The generalized Born rule asserts that if only the single Qbit is measured then the 1-Qbit measurement gate will indicate $`x`$ (0 or 1) with probability $`|a_x|^2`$, after which the $`n`$-Qbit state will be the product state $`|x|\mathrm{\Phi }_x`$, as illustrated in Figure 10.

Figure 10. Action of a 1-Qbit measurement gate on a single one of $`n`$ Qbits, according to the generalized Born rule.

To see that the gate acting on the measured Qbit in Figure 10 is indeed the $`n=1`$ version of the $`n`$-Qbit measurement gate of Figure 9, note that Figure 10 simplifies to Figure 11 in the special case $`|\mathrm{\Phi }_0=|\mathrm{\Phi }_1=|\mathrm{\Phi }`$. Now the entangled input state in Figure 10 becomes an uncorrelated product in which both the measured Qbit and the remaining $`n1`$ Qbits have states of their own. The $`(n1)`$ unmeasured Qbits now take no part whatever in the process. Nothing acts on them and they do nothing but maintain their original state $`|\mathrm{\Phi }`$. Their presence is irrelevant to the upper part of the figure, which is nothing more than the $`n=1`$ version of Figure 10.

Figure 11. Simplification of Figure 10 when $`|\mathrm{\Phi }_0=|\mathrm{\Phi }_1=|\mathrm{\Phi }`$

It is an elementary consequence of the generalized Born rule that the $`n`$-Qbit measurement gate of the ordinary Born rule can be constructed out of $`n`$ 1-Qbit measurement gates, as illustrated in Figure 12.

Figure 12. Constructing a 4-Qbit measurement gate out of four 1-Qbit measurement gates. The integer $`x`$ has the binary expansion $`x_3x_2x_1x_0`$.

Although the generalized Born rule is stronger, it follows from the Born rule under two plausible assumptions. (a) Once a Qbit ceases to interact with others and ceases to be acted on by unitary gates it doesn’t matter when it is measured. (b) To assign a state to Qbits is to do nothing more than to specify the probabilities of subsequent measurement outcomes, possibly preceded by unitary gates. Since the generalized Born rule reduces the notion of measurement to a single black-boxed 1-Qbit measurement gate (and indeed, since Qbits can be measured one by one, one needs only a single specimen of such a measurement gate) by far the most economical introduction to quantum computation is to base it on a primitive concept of the 1-Qbit measurement gate and make explicit assumptions (a) and (b) above.

Since it is impossible to determine the state of a collection of Qbits from the Qbits themselves, how is one to associate with the registers of a quantum computer the initial states on which the unitary transformations subsequently act? The obvious, simplest, and conceptually most economical answer is to exploit the measurement gates themselves. One can initialize $`n`$ Qbits to the state $`|\mathrm{\hspace{0.17em}0}_n`$ by measuring each Qbit and applying X if and only if the measurement indicates 1, as illustrated in Figure 13.

Figure 13. How to use a 1-Qbit measurement gate to prepare an off-the-shelf Qbit so its associated state is $`|\mathrm{\hspace{0.17em}0}`$.

Note the following features of this pedagogically motivated approach to quantum mechanics:

(1) It relies on an irreducible primitive notion of a unique black-boxed 1-Qbit measurement gate. The measurement gate is the only irreversible circuit element. It is defined by what it does, and what it does is to extract information from the Qbits in a form that is immediately accessible to us. There is no other way for us to get such information from the Qbits .

(2) Measurement plays a dual role. Output on a readable (classical) display not only ends the computation. It also provides, without any further complication, a straightforward way to begin the computation. Initialize every Qbit to the 1-Qbit state $`|\mathrm{\hspace{0.17em}0}`$ by sending each through a measurement gate, and then do nothing or apply X, depending on whether the display shows 0 or 1.

(3) Unlike the state of $`n`$ Cbits, the state of $`n`$ Qbits does not reside in the Qbits themselves: presented with a bunch of Qbits there is nothing one can do to them to reveal their state. Indeed, in general — for example if they share with other Qbits an entangled state — Qbits will not have a state of their own at all. To determine the state of Qbits (or whether they have one) one must go ask Alice, who knows their history: what initial measurements were performed on them, what the outcomes of the initial measurements were, and what subsequent unitary gates were applied.

(4) While the purpose of the state of $`n`$ Cbits is an anthropocentric add-on to its intrinsic character (what, for example, is the purpose of the velocity of a classical particle?) one would not bother to keep track of the state of $`n`$ Qbits were it not that this information about their past history has a specific purpose: it enables us to determine the correlations between the initial and final measurement outcomes after any intermediate sequence of applications of unitary gates.

All of these features resonate strongly with the constellation of ideas known as the Copenhagen interpretation. The quantum state of a system is not an objective property of that system. It merely provides an algorithm enabling one to infer from the initial set of measurements and outcomes (“state preparation”) the probabilities of the results of a final set of measurements after a specified intermediate time evolution. We ourselves have direct access to nothing beyond the outcomes of such measurements.

Why did my bare-bones, no-nonsense, pedagogically motivated, minimalist introduction to quantum mechanics come out sounding so Copenhagen? I think there are several reasons:

(a) Proponents of the Copenhagen interpretation (notably Heisenberg and Peierls) have always maintained that the quantum mechanical formalism does not describe “the system” but “our knowledge of the system”. Quantum computation is the first application of quantum mechanics that does not use it to further our understanding of the physical world. On the contrary, quantum computation exploits the known quantum mechanical character of the physical world to expedite the processing of knowledge, as represented symbolically by constituents (Qbits) of that world. It is therefore not surprising that the Copenhagen interpretation should provide a congenial setting for the exposition of quantum computation.

(b) A computation uses a finite set of Qbits. It has an unambiguous beginning and end. There is always a world external to the computation. If there were not an outside world there would be no point in doing the computation because there would be nobody or nothing to take advantage of the output. Nobody (well, at this stage practically nobody) wants to view the entire universe (single or multiple) as one colossal quantum computer, sufficient unto itself. The Copenhagen interpretation is characterized by a similar modesty of scope. Physics is a tool for relating some aspects of our experience to other aspects. Every application of physics begins and ends with an appeal to experience.

(c) The pedagogical device of restricting “measurement” to measurement in the computational basis, treating measurement in other bases as computational-basis measurement preceded by an appropriate unitary transformation, resonates with the Copenhagen notion of the primacy of the classical world. The computational basis states are (actually one should replace “are” by “are isomorphic to”) the states that describe ordinary classical Cbits. By restricting “measurement” to the computational basis I have automatically arranged for the input and output of every quantum computation to be describable in the ordinary old-fashioned language of classical computation — numbers flashed on a classical display. Bohr always insisted that our knowledge of the world had to be formulated in ordinary language, or we could not communicate it to anybody else.

(d) My computer science students know very little physics. They are therefore immune to any temptation to reify the states of Qbits into properties of the associated physical systems. If you think you too are immune from such temptation, ask yourself whether you do or do not belief that a horizontally polarized photon is intrinsically different from a vertically polarized photon. If you do, you are a victim of that very temptation.

I conclude by translating the possibly obscure quotation at the head of this section into the quite straightforward form it assumes in the context of quantum computation:

In our description of a quantum computation the purpose is not to disclose the real essence of the Qbits but only to track down statistical relations between initial and final measurement outcomes.

Acknowledgment

Supported by the National Science Foundation under grant No. 0098429.

References

1. D. DiVincenzo, P. Hayden, and B. Terhal, “Hiding quantum data”, www.arxiv.org/ abs/quant-ph/0207147.

2. C. H. Bennett et al., “Teleporting an unknown quantum state via dual classical and Einstein-Podolsky-Rosen channels”, Phys. Rev. Lett. 70, 1895-99 (1993).

3. C. A. Fuchs, “Notes on a Paulian idea: foundational, historical, anecdotal and forward-looking thoughts on the quantum”, www.arxiv.org/abs/quant-ph/0105039.

4. C. H. Bennet and S. J. Wiesner, “Communication via one- and two-particle operators on Einstein-Podolsky-Rosen states”, Phys. Rev. Lett. 69, 2881-84 (1992).

5. N. D. Mermin, “Deconstructing dense coding”, Phys. Rev. A 66, 032308 (2002).

6. I use here the unorthodox spelling Qbit because it will be constantly juxtaposed to Cbit, a role the currently favored qubit cannot gracefully play.

7. N. Bohr, Collected Works vol. 6, p. 296, North Holland (1985).

8. The pedagogical approach to quantum mechanics summarized here is described in more detail in N. D. Mermin, “From Cbits to Qbits: teaching computer scientists quantum mechanics”, American Journal of Physics 71, 23-30 (2003). Lecture notes for the course itself are regularly updated, revised, and reposted at http://www.ccmr.cornell .edu/~mermin/qcomp/CS483.html.

9. The term “c-bit” doesn’t work because one often needs to talk, for example, about 2-Cbit states.

10. The art of quantum computation, of course, is to construct unitary transformations leading to final states in which only informative values of $`x`$ are associated with appreciable probabilities $`|a_x|^2`$.


format: text

text:
# Glauber Dynamics on the Cycle is Monotone

## 1 Introduction

Consider the ferromagnetic Ising model on a finite graph $`G=(V,E)`$. Let $`n`$ denote the number of vertices. We allow variable interaction strength, so each edge $`e=[x,y]E`$ has a coupling constant $`J_{xy}>0`$. As usual, every vertex $`xV`$ has a spin $`\sigma _x\{+1,1\}`$ and the probability of a configuration is given by the Gibbs equilibrium measure $`\pi (\sigma )=Z^1\mathrm{exp}(_{[x,y]E}J_{xy}\sigma _x\sigma _y)`$. We assume no external field and no boundary conditions.

We are interested in the Glauber dynamics for this model. The Glauber dynamics is a Markov chain on the space of configurations of the Ising model, whose stationary distribution is the Gibbs measure. It comes in several flavors (heat-bath, Metropolis, etc.), but they follow the same scheme. At each step of the algorithm, one vertex is picked according to some rule (e.g. uniformly random) and with a certain probability one flips its spin (from $`+1`$ to $`1`$ and viceversa). The flip probability depends on the value of the spin and its neighbors; if it is chosen appropiately, then the chain distribution converges to the Gibbs measure.

For any configuration $`\sigma `$ we let $`\sigma ^x`$ be the configuration obtained by flipping spin $`x`$, while leaving all other spins unchanged. We shall focus on the heat-bath Glauber dynamics, in which a vertex $`x`$ is picked uniformly at random, then flipped with probability $`\pi (\sigma ^x)/(\pi (\sigma )+\pi (\sigma ^x))`$. Formally, this is the discrete-time Markov chain whose transition matrix $`A`$ has size $`2^n\times 2^n`$ and is defined by:

$`A(\sigma ,\sigma ^x)`$ $`=`$ $`(1/n)\pi (\sigma ^x)/(\pi (\sigma )+\pi (\sigma ^x))`$
$`A(\sigma ,\sigma )`$ $`=`$ $`1{\displaystyle \underset{xV}{}}A(\sigma ,\sigma ^x)`$
$`A(\sigma ,\tau )`$ $`=`$ $`0forallother\tau `$ (1)
One fundamental quantity is the relaxation time $`\tau _2`$, defined by

$$\tau _2=(n(1\mu _2))^1$$
(2)
where $`\mu _2`$ is the second largest eigenvalue of $`A`$. We scale by a factor of $`n`$ so that the $`\tau _2`$ defined above is the same as for a continuous-time chain where each vertex is updated at rate one. Essentially, it says we do $`n`$ moves of the discrete chain per unit of time.

In this paper we study the case when the graph $`G`$ consists of only one cycle (every vertex has degree two). Our main result is

###### Theorem 1

For the ferromagnetic Ising model on the cycle, the relaxation time $`\tau _2`$ is an increasing function of any of the couplings $`J_{xy}`$.

This result was motivated by trying to understand how the dynamics depends on the parameters of the model. The Glauber dynamics has been used in computer simulations of the Ising model, and in describing the non-equilibrium behavior of the model. There has been a fair amount of research on its relaxation time and mixing time; see  for a comprehensive survey, and  for an analysis of the Glauber dynamics on trees. Most existing results give asymptotic estimates for the relaxation time when the graph size is large; in particular, there are fairly sharp estimates for trees and for the lattice $`𝐙^d`$. However, much less is known about how the relaxation time depends on the parameters of the model, such as the temperature or the coupling constants. Such results would allow comparing mixing times for different graphs at different temperatures.

Computer simulations and heuristics suggest that, at least in the ferromagnetic case, such dependence should be monotone. Indeed, consider the dynamics started at an initial state where all spins are +1. To achieve mixing, lots of spins would have to flip to $`1`$. But in the ferromagnetic model, neighboring spins tend to have the same value. The higher the interaction strength, the harder it will be for spins to flip. This led to the following

###### Conjecture 2

(Yuval Peres, personal communication.) For the ferromagnetic Ising model on any graph, the relaxation time $`\tau _2`$ is an increasing function of any of the couplings $`J_{xy}`$.

What this says is that for two sets of couplings $`𝐉,𝐉^{}`$ with $`J_{xy}J_{xy}^{}`$, we must have $`\tau _2(𝐉)\tau _2(𝐉^{})`$. At first glance it may seem that one way to prove this is to find a coupling (not to be confused with the couplings $`J_{xy}`$) between the two Glauber dynamics, in which the dynamics for $`𝐉`$ is dominated by the one for $`𝐉^{}`$. This would imply the desired inequality for relaxation times.

However, such a coupling does not exist. To see why, choose some starting state (for example, all spins are +1) and let $`P_n,P_n^{}`$ be the measures obtained after $`n`$ steps of the dynamics. Let $`\pi ,\pi ^{}`$ be the corresponding Gibbs measures. We would like to define both dynamics on the same probability space so that stochastic dominance $`P_nP_n^{}`$ holds. But then it must also hold in the limit, so we would need to have $`\pi \pi ^{}`$, which cannot hold, because both measures have mean zero.

In this paper, we prove the conjecture in the special case when $`G`$ is a cycle. The proof proceeds as follows. We show that on the cycle, the transition matrix of the Glauber dynamics preserves linear functions of the spins, and furthermore that there is a linear eigenfunction corresponding to the second largest eigenvalue. Hence by looking at the action of the dynamics on the space of linear functions of spins, the original $`2^n\times 2^n`$ transition matrix can be reduced to a smaller, $`n\times n`$ matrix, whose entries are easily computable, and mostly zeroes.

We analyze the smaller matrix to obtain a variational description of the relaxation time, as well as an inequality involving the components of its lead eigenfunction. These two combine to prove the desired monotonicity. In addition, we obtain some further inequalities involving the relaxation time, as well as low-temperature asymptotics.

## 2 Two Useful Lemmas

The results in this section are valid for the ferromagnetic Ising model on any graph. They will be needed later in the proof of our main theorem.

The transition matrix $`A`$ of the heat-bath Glauber dynamics was defined above. It acts on the space of functions of $`\sigma `$ in the natural way. It has $`2^n`$ eigenvalues $`1=\mu _1>\mu _2\mathrm{}\mu _{2^n}`$ with corresponding eigenfunctions $`f_11,f_2,\mathrm{},f_{2^n}`$. Because $`A`$ is $`\pi `$-reversible we have $`\pi (f_if_j)=0`$ for $`ij`$ (if an eigenvalue is multiple, we can choose its eigenfunctions to be $`\pi `$-orthogonal), and in particular $`\pi f_i=0`$ for $`i2`$ (we use the standard notation $`\pi fE_\pi f`$).

An important class of functions are the increasing functions, defined in the usual way: $`f`$ is increasing if $`f(\sigma )f(\tau )`$ whenever $`\sigma \tau `$ (that is $`\sigma _i\tau _ii`$). $`f`$ is strictly increasing if strict inequality holds. We have

###### Lemma 3

(Y. Peres) The second eigenvalue $`\mu _2`$ has an increasing eigenfunction.

Proof. The method is similar to the proof of Perron-Frobenius. We need a standard result about the Glauber dynamics: if $`f`$ is increasing then $`Af`$ is increasing (for a proof see, for example, Theorem 2.14 in ).

Let $`f`$ be any increasing function with $`\pi f=0`$; it can be written as a linear combination of eigenvectors $`f=_{i=2}^{2^n}q_if_i`$. Let $`\mu =\mathrm{max}(|\mu _2|,|\mu _{2^n}|)`$. Iterating, we have

$$A^mf/\mu ^m=\underset{i=2}{\overset{2^n}{}}q_i(\mu _i/\mu )^mf_i$$
(3)
so if we let

$`g`$ $`=`$ $`{\displaystyle \underset{i=2}{\overset{2^n}{}}}q_if_i1(\mu _i=\mu )`$
$`h`$ $`=`$ $`{\displaystyle \underset{i=2}{\overset{2^n}{}}}q_if_i1(\mu _i=\mu )`$
we obtain convergence $`A^{2m}f/\mu ^{2m}g+h,A^{2m+1}f/\mu ^{2m+1}gh`$, where $`g`$ is some eigenvector corresponding to $`\mu _2`$ (a multiple of $`f_2`$ if $`\mu _2`$ is a simple eigenvalue). Since $`A^nf`$ is increasing for all $`n`$, the limits must also be increasing, so their sum must be increasing, so $`g`$ must be increasing.
It remains to show that $`g`$ is nonzero. If it were then we would have both $`h`$ and $`h`$ increasing, so $`h`$ would also have to be zero. To avoid this, it suffices to choose $`f`$ not orthogonal to $`f_2`$ and $`f_{2^n}`$. For example, $`f=f_2+f_{2^n}+K\sigma _i`$ has $`\pi f=0`$, and for large $`K>0`$, $`\pi (ff_2)0,\pi (ff_{2^n})0`$, and $`f`$ is increasing. Hence $`q_20,q_{2^n}0`$ so at least one of $`g`$ and $`h`$ is nonzero.

###### Lemma 4

(Y. Peres) If $`A`$ has a strictly increasing eigenfunction $`f`$, then $`f`$ corresponds to $`\mu _2`$.

Proof. Let $`g`$ be any other increasing eigenfunction. We will show there are only two possibilities: $`g`$ is constant (so its eigenvalue is 1), or $`g`$ corresponds to the same eigenvalue as $`f`$. Otherwise, we would have $`\pi (fg)=0,\pi f=0,\pi g=0`$. Since $`f`$ is strictly increasing, so is $`f\alpha g`$ for some $`\alpha >0`$ small enough. Since the Ising model has the FKG property (see ), $`\pi ((f\alpha g)g)\pi (f\alpha g)\pi g`$ so $`\pi (fg)\alpha \pi (g^2)(\pi f)(\pi g)\alpha (\pi g)^2`$ so $`\alpha \pi g^20`$, contradiction.
Lemma 3 guarantees there is an increasing $`g`$ corresponding to $`\mu _2`$; hence $`f`$ must also correspond to $`\mu _2`$.

## 3 Proof of the Main Result

From now on we assume our graph $`G`$ is a cycle, with vertices $`V=\{1,2,\mathrm{},n\}`$ and edges $`[1,2],[2,3],\mathrm{},`$ $`[n,1]`$. All notation will be mod $`n`$, so $`n+11`$. The coupling constants are $`J_{i,i+1}J_i0`$. As usual, every vertex $`iV`$ has a spin $`\sigma _i\{+1,1\}`$ and the probability of a configuration is given by $`\pi (\sigma )=Z^1\mathrm{exp}(_{i=1}^nJ_i\sigma _i\sigma _{i+1})`$. We assume no external field and no boundary conditions.

For technical reasons we assume $`J_i>0`$; the case when some of the $`J`$’s are zero can be easily obtained by limiting arguments. In particular, most results in this paper hold in the case when the graph $`G`$ is a straight line ($`J_{n,1}=0`$).

We have already introduced increasing functions. Another important class of functions are linear functions, that is functions of the form $`f(\sigma )=_{i=1}^nq_i\sigma _i`$ for some real constants $`q_i`$. On the cycle we have

###### Lemma 5

If $`f`$ is linear then $`Af`$ is linear.

Proof. We will need the standard hyperbolic functions $`\mathrm{sinh}(x)=1/2(\mathrm{exp}(x)\mathrm{exp}(x)),\mathrm{cosh}(x)=1/2(\mathrm{exp}(x)+\mathrm{exp}(x)),\mathrm{tanh}(x)=\mathrm{sinh}(x)/\mathrm{cosh}(x)`$. Since $`A`$ is a linear operator, it is enough to consider $`f(\sigma )=\sigma _i`$. Then $`f(\sigma ^j)=\sigma _i`$ for $`ji`$, and $`f(\sigma ^i)=\sigma _i`$, so

$`(Af)(\sigma )`$ $`=`$ $`A(\sigma ,\sigma )f(\sigma )+{\displaystyle \underset{j}{}}A(\sigma ,\sigma ^j)f(\sigma ^j)`$
$`=`$ $`(1{\displaystyle \underset{j}{}}A(\sigma ,\sigma ^j))\sigma _i+{\displaystyle \underset{j}{}}A(\sigma ,\sigma ^j)\sigma _i2A(\sigma ,\sigma _i)\sigma _i`$
$`=`$ $`\sigma _i2A(\sigma ,\sigma _i)\sigma _i`$
$`=`$ $`(11/n)\sigma _i+(1/n)\sigma _i(1/n){\displaystyle \frac{2\pi (\sigma ^i)}{\pi (\sigma )+\pi (\sigma ^i)}}\sigma _i`$
$`=`$ $`(11/n)\sigma _i+(1/n)\sigma _i{\displaystyle \frac{\pi (\sigma ^i)\pi (\sigma )}{\pi (\sigma ^i)+\pi (\sigma )}}`$
$`=`$ $`(11/n)\sigma _i+(1/n)\sigma _i\mathrm{tanh}(J_{i1}\sigma _i\sigma _{i1}+J_i\sigma _i\sigma _{i+1})`$
Since $`\mathrm{tanh}`$ is an odd function, the second term does not depend on $`\sigma _i`$, so

$$(Af)(\sigma )=(11/n)\sigma _i+(1/n)\mathrm{tanh}(J_{i1}\sigma _{i1}+J_{i+1}\sigma _{i+1})$$
(4)
and the $`\mathrm{tanh}`$ is easily expressed as a linear function

$$\mathrm{tanh}(J_{i1}\sigma _{i1}+J_{i+1}\sigma _{i+1})=a\sigma _{i1}+b\sigma _{i+1}$$
(5)
since $`\sigma _{i1}`$ and $`\sigma _{i+1}`$ only take the values $`+1`$ and $`1`$, and both sides of the equality are odd functions. So we only need to check that $`a+b=\mathrm{tanh}(J_{i1}+J_i),ab=\mathrm{tanh}(J_{i1}J_i)`$. After a bit of algebra, this yields

$`a`$ $`=`$ $`\mathrm{sinh}(2J_{i1})/(\mathrm{cosh}(2J_{i1})+\mathrm{cosh}(2J_i))`$
$`b`$ $`=`$ $`\mathrm{sinh}(2J_i)/(\mathrm{cosh}(2J_{i1})+\mathrm{cosh}(2J_i))`$
and the lemma is proven.

Hence we can look at the restriction of A to the space of linear functions. With respect to the basis $`\{\sigma _1,\sigma _2,\mathrm{},\sigma _n\}`$, this is described by an $`n\times n`$ matrix $`L`$ whose entries are all zero, except on the diagonal, and immediately above and below it. Let $`s_i=\mathrm{sinh}(2J_i)`$ and $`c_i=\mathrm{cosh}(2J_i)`$. Then we have

$`L(i,i)`$ $`=`$ $`11/n`$
$`L(i,i1)`$ $`=`$ $`(1/n)s_{i1}/(c_{i1}+c_i)`$
$`L(i,i+1)`$ $`=`$ $`(1/n)s_i/(c_{i1}+c_i)`$
Note that an eigenvalue of $`L`$ is also an eigenvalue of $`A`$, so all must be real. Since all entries on the diagonal of $`L`$ are equal, we can write $`L=(11/n)I+(1/n)M`$, where the diagonal of $`M`$ is all zeroes, and the only non-zero entries are

$$M(i,i1)=s_{i1}/(c_{i1}+c_i),M(i,i+1)=s_i/(c_{i1}+c_i)$$
(6)
Let the eigenvalues of $`M`$ be $`\lambda _1\mathrm{}\lambda _n`$. Then the eigenvalues of $`L`$ are $`1(1\lambda _i)/n`$ for $`1in`$, and $`M`$ and $`L`$ have the same eigenvectors. It is easy to see $`M`$ is irreducible. Hence by Perron-Frobenius, $`\lambda _1>0`$, $`\lambda _1`$ is a simple eigenvalue, and there exists a dominant eigenvector $`𝐱`$ (that is, corresponding to $`\lambda _1`$) which has non-negative entries. In fact, the entries must be strictly positive, since we have $`M𝐱=\lambda _1𝐱`$, so

$$M(i,i1)x_{i1}+M(i,i+1)x_{i+1}=\lambda _1x_i$$
(7)
so if $`x_i=0`$ for some $`i`$, then also $`x_{i+1}=0`$, so $`𝐱=0`$, contradiction. But then $`f(\sigma )=x_i\sigma _i`$ is a strictly increasing function of $`\sigma `$, and also an eigenfunction of $`A`$. Lemma 4 then implies that $`f`$ corresponds to the second eigenvalue. We have obtained

###### Corollary 6

$`\mu _2=1(1\lambda _1)/n`$, and $`\mu _2`$ has a linear eigenfunction.

Hence the relaxation time of the Glauber dynamics can be computed as $`\tau _2=1/(n(1\mu _2))=1/(1\lambda _1)`$, so from now on we can focus on the small matrix $`M`$. We will need an inequality which roughly says that its lead eigenvectors are not too far from a constant vector.

###### Lemma 7

If $`𝐱`$ is a non-negative dominant eigenvector of $`M`$, then

$$(c_i1)/s_ix_i/x_{i+1}(c_i+1)/s_i$$
(8)
Proof. The method could be called “cyclic induction.” Note that $`c_i^2s_i^2=1`$, so $`(c_i1)/s_i=s_i/(c_i+1)<1`$. We have $`M𝐱=\lambda _1𝐱`$, so

$$\lambda _1(c_{i1}+c_i)x_i=s_{i1}x_{i1}+s_ix_{i+1}i$$
(9)
Assume the left-hand inequality of (8) holds for $`i1`$ so $`(c_{i1}1)/s_{i1}x_{i1}/x_i`$. Combining this with (9), we get (since $`\lambda _1<1`$)

$$(c_{i1}+c_i)x_i\lambda _1(c_{i1}+c_i)x_i(c_{i1}1)x_i+s_ix_{i+1}$$

so $`(c_i+1)x_is_ix_{i+1}`$ so $`x_i/x_{i+1}s_i/(c_i+1)=(c_i1)/s_i`$. Hence the inequality holds for $`i`$ as well; hence it holds either for all $`i`$ or for none. But in the latter case we would get $`x_i/x_{i+1}<(c_i1)/s_i<1`$ for all $`i`$, hence $`x_1<x_2<\mathrm{}<x_n<x_1`$, which is impossible.

The right-hand inequality can be proven in a similar way: if the inequality $`(c_i+1)/s_ix_i/x_{i+1}`$ holds for $`i1`$, it also holds for $`i`$; it cannot hold for all $`i`$, so it cannot hold for any $`i`$. Alternatively, the right-hand side can be deduced from the left-hand side by symmetry.

An immediate consequence of this is

###### Proposition 8

If $`𝐱`$ is a non-negative dominant eigenvector of $`M`$, then $`x_i/x_{i+1}+x_{i+1}/x_i2c_i/s_i`$.

Proof. The function $`f(t)=t+1/t`$ is decreasing for $`t<1`$ and increasing for $`t>1`$, so it follows from Lemma 7 that $`f(x_i/x_{i+1})f((c_i+1)/s_i)`$ if $`x_ix_{i+1}`$, and $`f(x_i/x_{i+1})f((c_i1)/s_i)`$ if $`x_ix_{i+1}`$. Both upper bounds are equal to $`2c_i/s_i`$.

To apply the last result we need an expression for the largest eigenvalue in terms of its eigenvector. Note that the product of the elements under the diagonal of $`M`$ is equal to the product of the ones above the diagonal, so $`M`$ is reversible in the following sense. Define $`q(i)=c_i+c_{i+1}`$. Then $`q(i)M(i,j)=q(j)M(j,i)i,j`$. So if $`Q`$ is the diagonal matrix with $`Q(i,i)=q(i)`$ then $`Q^{1/2}MQ^{1/2}`$ is symmetric, and has the same eigenvalues as $`M`$. Then by the usual extremal characterization,

$$\lambda _1=sup\frac{(Q^{1/2}MQ^{1/2}v,v)}{(v,v)}$$

where the $`sup`$ can be taken over all nonzero vectors $`v`$, and is attained at dominant eigenvectors of $`Q^{1/2}MQ^{1/2}`$. The numerator expands as

$`(Q^{1/2}MQ^{1/2}v,v)`$ $`=`$ $`{\displaystyle 2v_iv_{i+1}\sqrt{M(i,i+1)M(i+1,i)}}`$
$`=`$ $`{\displaystyle 2v_iv_{i+1}s_i/\sqrt{(c_{i1}+c_i)(c_i+c_{i+1})}}`$
so if we finally substitute $`v_i=(c_i+c_{i+1})^{1/2}y_i`$ we obtain

###### Lemma 9

$`\lambda _1=sup_𝐲(2y_iy_{i+1}s_i)/((y_i^2+y_{i+1}^2)c_i)`$ and the $`sup`$ is attained at dominant eigenvectors of $`M`$.

Now let $`h(𝐲,𝐉)=(2y_iy_{i+1}s_i)/((y_i^2+y_{i+1}^2)c_i)`$ and let $`𝐱𝐱(𝐉)`$ be a non-negative dominant eigenvector of $`M`$, chosen so $`𝐱()`$ is a smooth function in a neighborhood of $`𝐉`$. This can be done because $`\lambda _1`$ is a simple eigenvalue; see, for example, , Chapter 2.1. We have $`\lambda _1(𝐉)=h(𝐱(𝐉),𝐉)`$. Because for fixed $`𝐉`$, $`suph(𝐲,𝐉)`$ is attained at $`𝐲=𝐱(𝐉)`$, we have

$$h/y_k(𝐱(𝐉),𝐉)=0k$$
(10)
so by the chain rule

$$\lambda _1/J_i=h/J_i+\underset{k}{}(h/y_k)(x_k/J_i)=h/J_i$$
(11)
Hence

$`\lambda _1/J_i`$ $`=`$ $`h/J_i(𝐱,𝐉)`$ (12)
$`=`$ $`2(2x_ix_{i+1}c_i(x_i^2+x_{i+1}^2)s_ih(𝐱,𝐉))/{\displaystyle (x_i^2+x_{i+1}^2)c_i}`$
From Proposition 8, $`(x_i^2+x_{i+1}^2)/(2x_ix_{i+1})c_i/s_i`$ and since $`h(𝐱,𝐉)=\lambda _1<1`$ it follows that $`\lambda _1/J_i>0`$. This completes the proof of Theorem 1.

## 4 Further Inequalities

We can get more precise bounds on the derivatives of $`\lambda _1`$. Indeed, equation (12) and Proposition 8 yield easily

$$\lambda _1/J_i2(1\lambda _1)(2x_ix_{i+1}c_i)/(x_i^2+x_{i+1}^2)c_i$$
(13)
hence

$`{\displaystyle (s_i/c_i)\lambda _1/J_i}`$ $``$ $`2(1\lambda _1){\displaystyle (2x_ix_{i+1}s_i)/(x_i^2+x_{i+1}^2)c_i}`$ (14)
$`=`$ $`2(1\lambda _1)\lambda _1`$
So we can obtain an interesting inequality by differentiating along a particular direction. To do that, consider functions $`J_i(t)`$ with $`J_i(0)=J_i`$ and $`J_i^{}(t)=\mathrm{tanh}(2J_i(t))`$. This differential equation is easy to solve and yields

$$J_i(t)=(1/2)\mathrm{sinh}^1(\mathrm{exp}(2t)\mathrm{sinh}(2J_i))$$
(15)
(here $`\mathrm{sinh}^1`$ denotes the inverse hyperbolic function). Then the function $`\lambda (t)=\lambda _1(J_1(t),\mathrm{},J_n(t))`$ satisfies $`\lambda ^{}2\lambda (1\lambda )`$. This becomes even simpler if we substitute

$$\tau ^{}=\tau _21=\lambda /(1\lambda )$$
(16)
We obtain $`\tau _{}^{}{}_{}{}^{}2\tau ^{}`$ so $`(\mathrm{log}\tau ^{})^{}2`$ and $`\tau ^{}(t)\mathrm{exp}(2t)`$. Hence, with a slight abuse of notation (by $`f(𝐉)`$ we mean that $`f`$ is applied to all entries of $`𝐉`$), we have

###### Proposition 10

$`\tau ^{}((1/2)\mathrm{sinh}^1(\mathrm{exp}(2t)\mathrm{sinh}(2𝐉)))\mathrm{exp}(2t)\tau ^{}(𝐉)𝐉,t0`$.

For large $`x`$, $`\mathrm{sinh}^1(x)`$ is close to $`\mathrm{log}(x)`$, so we expect $`J_i(t)`$ to be approximately $`J_i+t`$. In fact it is easy to prove that $`\mathrm{exp}(2t)\mathrm{sinh}(2u)\mathrm{sinh}(2t+2u)`$, so $`J_i(t)J_i+t`$, so from Theorem 1 we get

###### Theorem 11

$`\tau ^{}(𝐉+t)\mathrm{exp}(2t)\tau ^{}(𝐉)`$ for all $`t0`$.

As above, $`𝐉+t`$ means $`t`$ is added to all entries of $`𝐉`$. Numerical computations confirm that the factor 2 in the exponent is optimal. It would be nice to have a similar inequality where the elements of $`𝐉`$ are multiplied by $`t`$, so that we can estimate the effect of a change of temperature on the relaxation time. Note, however, that $`t𝐉>t+𝐉`$ for large $`t`$, so the inequality we have proven is in some sense stronger.

## 5 Low Temperature Asymptotics

Introduce the inverse temperature $`\beta `$; the coupling constants are now $`\beta J_1,\mathrm{},\beta J_n`$. We are interested in $`\tau _2(\beta 𝐉)`$ for large $`\beta `$ and fixed $`𝐉`$. This turns out to be asymptotically exponential in $`\beta `$. More precisely, let $`J_{max}=\mathrm{max}J_i`$ and $`J_{min}=\mathrm{min}J_i`$. We have

###### Proposition 12

The following holds as $`\beta \mathrm{}`$:

$$lim\tau _2(\beta 𝐉)/\mathrm{exp}(C\beta )=D$$

where $`C=2(J_{max}+J_{min})`$ and $`D=(1/2)\mathrm{\#}\{i:J_i=J_{max}\}/\mathrm{\#}\{i:J_i=J_{min}\}`$.

Proof. As before, let $`𝐱`$ be a dominant eigenvector. Lemma 7 guarantees that the ratios $`x_i/x_{i+1}`$ are close to 1 for large $`\beta `$. That is, given $`ϵ>0`$, $`1ϵx_i/x_{i+1}1+ϵ`$ for $`\beta `$ large enough. We can set $`x_1=1`$ so we get $`(1ϵ)^nx_i(1+ϵ)^ni`$.

If we apply Lemma 9 for $`𝐲1`$ we obtain $`\lambda _1s_i/c_i`$ so

$$\tau _2=1/(1\lambda _1)c_i/(c_is_i)$$
(17)
If we apply Lemma 9 for $`𝐲=𝐱`$ and substitute $`\tau _2=1/(1\lambda _1)`$, we obtain

$`\tau _2`$ $`=`$ $`{\displaystyle c_i(x_i^2+x_{i+1}^2)/[c_i(x_ix_{i+1})^2+2x_ix_{i+1}(c_is_i)]}`$ (18)
$``$ $`((1+ϵ)/(1ϵ))^{2n}{\displaystyle c_i/(c_is_i)}`$
The ratio $`c_i/(c_is_i)`$ is easily shown to be asymptotically $`D\mathrm{exp}(C\beta )`$. Since $`ϵ`$ is arbitrary, the proof is complete.

Acknowledgments. Many thanks to Yuval Peres, my thesis advisor. His guidance was instrumental in obtaining these results. I also thank him for permission to include Lemmas 3 and 4. Thanks also to Noam Berger, Elchanan Mossel and Dror Weitz for their comments.


text:
# Streamer Branching rationalized by Conformal Mapping Techniques

## Abstract

Spontaneous branching of discharge channels is frequently observed, but not well understood. We recently proposed a new branching mechanism based on simulations of a simple continuous discharge model in high fields. We here present analytical results for such streamers in the Lozansky-Firsov limit where they can be modelled as moving equipotential ionization fronts. This model can be analyzed by conformal mapping techniques which allow the reduction of the dynamical problem to finite sets of nonlinear ordinary differential equations. Our solutions illustrate that branching is generic for the intricate head dynamics of streamers in the Lozansky-Firsov-limit.

When non-ionized matter is suddenly exposed to strong fields, ionized regions can grow in the form of streamers. These are ionized and electrically screened channels with rapidly propagating tips. The tip region is a very active impact ionization region due to the self-generated local field enhancement. Streamers appear in early stages of atmospheric discharges like sparks or sprite discharges , they also play a prominent role in numerous technical processes. It is commonly observed that streamers branch spontaneously . But how this branching is precisely determined by the underlying discharge physics, is essentially not known. In recent work , we have suggested a branching mechanism from first principles. This work drew some attention , since the proposed mechanism yields quantitative predictions for specific parameters, and since it is qualitatively different from the older branching concept of the “dielectric breakdown model” . This older concept actually can be traced back to concepts of rare long-ranged (and hence stochastic) photo-ionization events probably first suggested in 1939 by Raether . Therefore, it came as a surprise that we predicted streamer branching in a fully deterministic model. Since our evidence for the phenomenon was mainly from numerical solutions together with a physical interpretation, the accuracy of our numerical scheme was challenged . Furthermore, some authors have argued previously  that in a deterministic discharge model like ours, an initially convex streamer head never could become locally concave, and that hence the consecutive branching of the discharge channel would be unphysical.

Therefore in the present paper, we investigate the issue by analytical means. We show that the convex-to-concave evolution of the streamer head with successive branching is generic for streamers in the Lozansky-Firsov limit . We define the Lozansky-Firsov-limit as the stage of evolution where the streamer head is almost equipotential and surrounded by a thin electrostatic screening layer. While in the original article , only simple steady state solutions with parabolic head shape are discussed, we will show here that a streamer in the Lozansky-Firsov limit actually can exhibit a very rich head dynamics that includes spontaneous branching. Furthermore, our analytical solutions disprove the reasoning of  by explicit counterexamples. Our analytical methods are adapted from two fluid flow in Hele-Shaw cells . But our explicit and exact solutions that amount to the evolution of “bubbles” in a dipole field , have not been reported in the hydrodynamics literature either.

The relation between our previous numerical investigations  and our present analytical model is laid in two steps. First, numerical solutions show essentially the same evolution in the purely two-dimensional case as in the three-dimensional case with assumed cylinder geometry . Because there is an elegant analytical approach, we focus on the two-dimensional case. This has the additional advantage that such two-dimensional solutions rather directly apply to, e.g., discharges in Corbino discs . Second, we use the following simplifying approximations for a Lozansky-Firsov streamer: $`(1)`$ the interior of the streamer is electrically completely screened; hence the electric potential $`\phi `$ is constant; hence the ionization front coincides with an equipotential line, $`(2)`$ the width of the screening layer around the ionized body is much smaller than all other relevant length scales and in the present study it is actually neglected, $`(3)`$ the velocity of the ionization front v is determined by the local electric field; in the simplest case to be investigated here, it is simply taken to be proportional to the field at the boundary $`𝐯=c\phi `$ with some constant $`c`$ (for the validity of the approximation, cf. ). Together with $`^2\phi =0`$ in the non-ionized outer region and with fixed limiting values of the potential $`\phi `$ far from the streamer, this defines a moving boundary problem for the interface between ionized and non-ionized region. We assume the field far from the streamer to be constant as in our simulations . Such a constant far field can be mimicked by placing the streamer between the two poles of an electric dipole where the distance between the poles is much larger than the size of the streamer.

When the electric field points into the $`x`$ direction and $`y`$ parametrizes the transversal direction, our two-dimensional Lozansky-Firsov streamer in free flight in a homogeneous electric field is approximated by:

$`^2\phi (x,y)=0,`$ $`\text{outside the streamer},`$ (1)
$`\phi (x,y)E_0\widehat{x},`$ $`\text{far outside the streamer},`$ (2)
$`\phi (x,y)=0,`$ $`\text{inside the streamer},`$ (3)
$`𝐯_{\mathrm{bound}.}=c\phi |_{\mathrm{bound}.},`$ $`\text{velocity of the boundary},`$ (4)
where $`\widehat{x}`$ is the unit vector in the $`x`$ direction, and we have chosen the gauge such that the potential inside the streamer vanishes. The asymptote (2) implies that the total charge on the streamer vanishes; otherwise a contribution $`1/\sqrt{x^2+y^2}`$ has to be added on the r.h.s. of Eq. (2).

Similar moving boundary problems arise in Hele-Shaw flow of two fluids with a large viscosity contrast : Lozansky-Firsov streamers and viscous fingers on the present level of description can be identified immediately by equating the electric potential $`\phi `$ with the fluid pressure $`p`$ . To such problems, powerful conformal mapping methods  can be applied. Most work with this method is concerned with viscous fingers in a channel geometry, i.e., with boundary conditions on a lateral external boundary that cannot be realized in an electric system. A few authors also study air bubbles within a viscous fluid, or viscous droplets in air, mostly under the action of flow fields generated by one source or one sink of pressure, i.e., by monopoles. On the other hand, the approximation (1)-(4) describes streamers in free space between two electrodes as in . With the asymptote (2), this is mathematically equivalent to air bubbles in a dipole field. This case has not been studied in detail. It is known that any ellipse with the main axes oriented parallel and perpendicular to the direction of the dipole is a uniformly translating solution of this problem . The time dependent solutions of  do not apply to streamers since the boundary condition on the moving interface is different. Refs.  and  study how and when cusps in the interfaces of droplets and bubbles emerge when these are driven by multipole fields. But for bubbles in a dipole field, again only the steady state ellipse solutions are given .

In the present paper, we therefore apply conformal mapping methods to the evolution of “bubbles” in a dipole field in a Hele-Shaw experiment and proceed beyond the steady state ellipse solutions. We identify the general structure of time dependent solutions of (1)-(4). The analytically derived exact solutions show how a streamer head can become flatter, concave and branch as observed numerically . Rather than a pole decomposition , we derive a decomposition into Fourier modes of the circle and calculate an equation for the non-linear dynamical coupling of their amplitudes.

In detail, this is done in the following steps:
$`(i)`$ The spatial coordinates are expressed by the complex coordinate $`z=x+iy`$. According to standard complex analysis, finding a real harmonic function $`\phi (x,y)`$ solving the Laplace equation (1) in a given region is equivalent to finding a complex function $`\mathrm{\Phi }(z)`$ that is analytical in the same region and has real part $`\mathrm{Re}\mathrm{\Phi }(z)=\phi (x,y)`$.
$`(ii)`$ A conformal map from the interior of the unit circle to the exterior of the streamer or “bubble” is constructed. Including the point at infinity, the region outside the streamer is simply connected and Riemann’s mapping theorem applies; therefore the mapping exists. Since the boundary moves, the mapping is time dependent; we denote it with $`z=f_t(\omega )`$ where $`\omega `$ parametrizes the interior of the unit circle $`|\omega |<1`$. The complete map can be composed from a conformal map $`\zeta =h_t(\omega )`$ that deforms the unit disc continuously, followed by the inversion $`z=1/\zeta `$. Since $`h_t(\omega )`$ is conformal on the unit disc, it is analytical and has a single zero which we choose to be at $`\omega =0`$. Therefore $`f_t(\omega )=1/h_t(\omega )`$ has a single pole $`\omega ^1`$ and is otherwise analytical. Rather than a pole decomposition , we choose a Laurent expansion for $`f_t(\omega )`$

$`x+iy=z=f_t(\omega )={\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}a_k(t)\omega ^k.`$ (5)
This expansion allows us to identify the exact dynamical solutions (12) below. Taking $`a_1(t)`$ as a real positive number makes the mapping unique, again according to Riemann’s mapping theorem.
$`(iii)`$ Now the potential $`\widehat{\mathrm{\Phi }}(\omega )`$ on the unit disc can be calculated explicitly. Since $`f_t(\omega )`$ is a conformal mapping, the function $`\mathrm{\Phi }(z)`$ is analytical if and only if the function $`\widehat{\mathrm{\Phi }}(\omega )=\mathrm{\Phi }(f_t(\omega ))`$ is analytical. The asymptote of $`\widehat{\mathrm{\Phi }}(\omega )`$ for $`\omega 0`$ is determined by (2) and (5): for $`|x|,|z|\mathrm{}`$, we have $`\phi (x,y)E_0x`$, hence $`\mathrm{\Phi }(z)E_0z`$, and therefore with (5): $`\widehat{\mathrm{\Phi }}(\omega )E_0a_1(t)/\omega `$ for $`\omega 0`$. This means that the pole of $`\widehat{\mathrm{\Phi }}(\omega )`$ at the origin of the unit disc $`\omega =0`$ corresponds to the dipole of $`\mathrm{\Phi }(z)`$ at $`z\pm \mathrm{}`$. This dipole generates the field and the interfacial motion. In the remainder of the unit disc, there are no sources or sinks of potential, hence $`\widehat{\mathrm{\Phi }}`$ is analytical there. Furthermore, at the boundary of the streamer, we have $`\phi =0`$ from (3) or $`\mathrm{Re}\mathrm{\Phi }=0`$, resp. The boundary of the streamer maps onto the unit circle, so $`\mathrm{Re}\widehat{\mathrm{\Phi }}(\omega )=0`$ for $`|\omega |=1`$. Using the asymptotics at $`\omega 0`$ and analyticity in the remaining region, the unique and exact solution for the potential is

$`\widehat{\mathrm{\Phi }}(\omega )=E_0a_1(t)\left(\omega {\displaystyle \frac{1}{\omega }}\right).`$ (6)
$`(iv)`$ The velocity $`𝐯_{\mathrm{bound}.}=c\varphi `$ (4) determines the motion of the interface. This interface is the time dependent map $`f_t(\omega )`$ of the unit circle $`\omega =e^{i\alpha }`$ parametrized by the angle $`\alpha [0,2\pi )`$. Therefore Eq. (4) determines the dynamics of $`f_t\left(e^{i\alpha }\right)`$. According to , it is

$$\mathrm{Re}\left[i_\alpha f_t^{}\left(e^{i\alpha }\right)_tf_t\left(e^{i\alpha }\right)\right]=c\mathrm{Re}\left[i_\alpha \widehat{\mathrm{\Phi }}\left(e^{i\alpha }\right)\right].$$
(7)
The problem (1)-(4) is symmetric under reflection on the $`x`$-axis. We assume the solutions to have the same symmetry. Therefore all $`a_k(t)`$ have to be real. The position $`(x,y)(\alpha ,t)`$ of the point of the interface labelled by the angle $`\alpha `$ at time $`t`$ can be read directly from the Laurent expansion (5) by inserting $`\omega =e^{i\alpha }`$; it has essentially the form of a Fourier expansion of the unit circle where the circle and its position are created by the modes $`k=1`$ and 0:

$`x(\alpha ,t)={\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}a_k(t)\mathrm{cos}k\alpha ,y(\alpha ,t)={\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}a_k(t)\mathrm{sin}k\alpha ,`$ (8)
$`a_k(t)\text{real},a_1(t)>0,0\alpha <2\pi .`$ (9)
Substituting the mapping function (5) and the potential (6) into the equation of motion for the mapping (7), and assuming the $`a_k(t)`$ to be real, we obtain for the evolution of the amplitudes $`a_k(t)`$:

$`{\displaystyle \underset{k,k^{}=1}{\overset{\mathrm{}}{}}}k^{}a_k^{}(t)_ta_k(t)\mathrm{cos}\left((kk^{})\alpha \right)`$ (10)
$`=2E_0ca_1(t)\mathrm{cos}\alpha .`$ (11)
A closer investigation shows that this equation has an important property: suppose that the streamer boundary can be written initially as a finite series $`_{k=1}^Na_k(0)e^{ik\alpha }`$, $`a_N(0)0`$. Then at all times $`t`$, the interface is described by the same finite number of modes

$`z(\alpha ,t)={\displaystyle \underset{k=1}{\overset{N}{}}}a_k(t)e^{ik\alpha },`$ (12)
i.e., the $`a_k(t)`$ with $`k>N`$ stay identical to zero at all times $`t>0`$. Sorting the terms in (10) by coefficients of $`\mathrm{cos}k\alpha `$, the equation can be recast into $`N+2`$ ordinary differential equations for the $`N+2`$ functions $`a_k(t)`$

$`{\displaystyle \underset{k=1}{\overset{Nm}{}}}\left[(k+m)a_{k+m}_ta_k+ka_k_ta_{k+m}\right]`$ (14)
$`=2E_0ca_1\delta _{m,1}\text{ for }m=0,\mathrm{},N+1,`$
where $`\delta _{m,1}`$ is the Kronecker symbol. Eq. (14) is equivalent to a matrix equation of the form $`𝐀(\{a_k(t)\})_t(a_1(t),\mathrm{},a_N(t))=(0,2E_0ca_1(t),0,\mathrm{},0)`$, where the matrix $`𝐀`$ depends linearly on the $`\{a_k(t)\}`$.

Eqs. (12) and (14) identify large classes of analytical solutions with arbitrary fixed $`N`$. These solutions reduce the dynamical moving boundary problem in two spatial dimensions of Eqs. (1)–(4) exactly to a finite set of ordinary differential equations for the nonlinear coupling of the amplitudes $`a_k(t)`$ of modes $`e^{ik\alpha }`$, $`0\alpha <2\pi `$. These equations are easy to integrate numerically or for small $`N`$ even analytically. We will use this form to discuss now generic solutions of Eqs. (1)–(4) as the simplest approximation of a streamer in the Lozansky-Firsov limit.

First, it is now easy to reproduce the uniformly propagating ellipse solutions of  as the solutions with $`N=1`$: for $`|a_1||a_1|`$, the equations reduce to $`_ta_1=0=_ta_1`$ and $`_ta_0=2E_0ca_1/(a_1a_1)`$. These solutions correspond to ellipses whose principal radii are oriented along the axes. These radii maintain their values $`r_{x,y}=a_1\pm a_1`$ (assuming $`a_1>a_1>0`$) and move with constant velocity $`v_{\mathrm{ellipse}}=E_0c(r_x+r_y)/r_y`$. The Lozansky-Firsov-parabola can be understood as limit cases of such uniformly propagating ellipses.

In contrast to $`N1`$, all solutions with $`N2`$ have nontrivial dynamics. It can be tracked by integrating the $`N+2`$ ordinary differential equations (14) numerically and then plotting the boundaries (12) at consecutive times. Examples of such dynamics are shown in the figures.

Fig. 1 shows four cases of the upward motion of a conically shaped streamer in equal time steps. The initial conditions are almost identical. On the leftmost figure, an ellipse is corrected only by a mode $`e^{2i\alpha }`$ to create the conical shape. This shape with $`N=2`$ eventually develops a concave tip, but only after much longer times than shown in the figure. In the other figures this conical shape is perturbed initially by a minor perturbation with wavenumber 8 or 30, corresponding to $`N=8`$ and 30 in (12) and (14). The amplitude of the perturbation is chosen such that a cusp develops at time $`0.1/(E_0c)`$. Depending on the sign of the amplitude, the cusp develops on or off axis — where we stress that we are not interested in the cusp itself, but in the earlier stages of evolution. Note that our reduction of the moving boundary problem to the set of ordinary differential equations (14) assures that the evolving shape is a true solution of the problem (1)–(4). Figs. $`1b`$ and $`1d`$ demonstrate that spontaneous branching is a possible solution.

In Fig. 2 the ionized body is longer stretched and only the tip is shown, again at 6 equidistant time steps. The streamer becomes slower when the head becomes flatter, since the electric field then diminishes. Eventually, the head becomes concave and “branches”.

In summary, the solutions of the moving boundary problem (1)–(4) demonstrate the onset of branching within a purely deterministic model. They show a high sensitivity to minor deviations of the initial conditions. A streamer in the Lozansky-Firsov-limit is therefore also very sensitive to physical perturbations during the evolution, and simulations in this limit are just as sensitive to small numerical errors. But perturbations during the evolution are not necessary for branching.

Our analysis applies to streamers in the Lozansky-Firsov-limit, i.e., to almost equipotential streamers that are surrounded by a very thin electrical screening layer. This limit is approached in our previous simulations .

These results raise the following questions that are presently under investigation: 1) When does a streamer reach this Lozansky-Firsov-limit that then generically leads to branching? 2) The formation of cusps should be suppressed by some microscopic stabilization mechanism. Is the electric screening length discussed in  sufficient to supply this mechanism? 3) If this stabilization is taken into account, can an interfacial model reproduce numerical and physical streamer branching quantitavely? 4) How can the motion of the back end of the streamer be modelled appropriately (rather than assuming the velocity law $`𝐯\phi `$ (4) everywhere)? How can it be incorporated into the present analysis?

Acknowledgment: B.M. was supported by CWI Amsterdam, and A.R. by the Dutch Research School CPS, the Dutch Physics Funding Agency FOM and CWI.


text:
# Total Reaction Cross Section in an Isospin-Dependent Quantum Molecular Dynamics (IDQMD) Model

thanks:  Supported by the Major State Basic Research Development Program of China under Grant No TG 2000077404.

## Abstract

The isospin-dependent quantum molecular dynamics (IDQMD) model is used to study the total reaction cross section $`\sigma _R`$. The energy-dependent Pauli volumes of neutrons and protons have been discussed and introduced into the IDQMD calculation to replace the widely used energy-independent Pauli volumes. The modified IDQMD calculation can reproduce the experimental $`\sigma _R`$ well for both stable and exotic nuclei induced reactions. Comparisons of the calculated $`\sigma _R`$ induced by $`{}_{}{}^{11}Li`$ with different initial density distributions have been performed. It is shown that the calculation by using the experimentally deduced density distribution with a long tail can fit the experimental excitation function better than that by using the Skyrme-Hartree-Fock calculated density without long tails. It is also found that $`\sigma _R`$ at high energy is sensitive to the long tail of density distribution.

The total reaction cross section $`\sigma _R`$ has been extensively studied theoretically and experimentally.\[1-10\] There are two kinds of theoretical models to calculate $`\sigma _R`$. The first is the low energy theory based on the interaction potential. Such models are not successful for the reaction beyond 10-15 MeV/nucleon above the Coulomb barrier. The second is the high energy microscopic Glauber theory based on the individual nucleon-nucleon collisions in the overlap volume of the projectile and target. However, the roles of mean field and medium effect are difficult to be discussed in the Glauber-type model. In recent years, Ma et al. developed a new method to study $`\sigma _R`$ with helps of transport theory and Glauber model. Originally, the BUU model  was taken as a tool to investigate $`\sigma _R`$. Later, the quantum molecular dynamics (QMD) model was applied to study $`\sigma _R`$  in the same spirit as Ref. . The reaction dynamics in transport theory at intermediate energy is mainly governed by the mean field, two-body collisions, and Pauli blocking. To investigate the isospin effects, the above three dynamical ingredients should include properly isospin degrees of freedom to obtain an isospin-dependent quantum molecular dynamics (IDQMD). It is also important that the samples of neutrons and protons in the phase space should be treated separately in the initialization of projectile and target nuclei. In this Letter, the IDQMD model is introduced to calculate $`\sigma _R`$. This model incorporates the isospin-dependence of mean-field, nucleon-nucleon cross section, and Pauli blocking. It has been widely used to study the multi-fragmentation and the collective flow. In the IDQMD model, neutrons and protons are distinguishable in the initialization.

The density distributions for the initial projectile and target nuclei are determined from the Skyrme-Hartree-Fock (SHF) calculation  with parameter set SKM. The stability of the propagation of the initialized nuclei has been checked in details and can last at least 200fm/c according to the evolutions of the average binding energies and the root mean square radii of the initialized nuclei.

In the IDQMD model, once the distance between the two nucleons is less than rnn, the collision may occur, where rnn is defined as

$$r_{nn}=\sqrt{\sigma _{nn}(\sqrt{s})/\pi }$$

with $`\sigma _{nn}(\sqrt{s})`$ representing the nucleon-nucleon reaction cross section and $`\sqrt{s}`$ representing the nucleon energy in the center of mass system. Whenever a collision occurs, the final momenta of the scattering nucleons can be easily obtained as a consequence of momentum conservation and the coordinates are updated in terms of newtonian equation (classical trajectory). The six-dimensional phase space radius of one nucleon is the product of its $`\mathrm{\Delta }`$P and $`\mathrm{\Delta }`$R. Then we calculate and check the phase spaces around the final states of the scattering nucleons. Thus it is easy to determine the probabilities ($`P_1`$ and $`P_2`$) for each of the two scattering nucleons that their final phase spaces are already occupied by other nucleons by comparing with the defined Pauli volue, i.e. 4$`(\mathrm{\Delta }P\mathrm{\Delta }R)^3/3`$, in six-dimensional phase space, where $`\mathrm{\Delta }R`$ is the minimum radius which is allowed to be occupied by itself in corrdinate space and the $`\mathrm{\Delta }P`$ is the same but in momentum space. The collision is then blocked with a probability

$$P_{block}=1[1min(P_1,1)][1min(P_2,1)].$$

Correspondingly, the collision is allowed with probability $`(1P_{block})`$. Whenever a collision is blocked, the momenta of the scattering nucleons will be replaced by the values they are prior to scattering. More details could be found in Refs. .

The following formula of in-medium reaction cross section is used in the present IDQMD calculation,

$`\sigma _{nn}=(13.7315.04\beta ^1+8.76\beta ^2+68.67\beta ^4)`$
$`{\displaystyle \frac{1.0+7.772E_{lab}^{0.006}\rho ^{1.48}}{1.0+18.0\rho ^{1.46}}}`$
$`\sigma _{np}=(70.6718.18\beta ^1+25.26\beta ^2+113.85\beta )`$
$`{\displaystyle \frac{1.0+20.88E_{lab}^{0.04}\rho ^{2.02}}{1.0+35.86\rho ^{1.90}}}`$
$`\beta =\sqrt{1{\displaystyle \frac{1.0}{\gamma ^2}}},\gamma ={\displaystyle \frac{E_{lab}}{931.5}}+1`$ (1)
This formula includes both effects of incident energy Elab and nucleon matter density $`\rho `$. It shows that the medium effect is important at intermediate energies and becomes smaller at higher energies but does not vanish..

The cross section $`\sigma _R`$ can be written as

$$\sigma _R=2\pi b[1T(b)]𝑑b=2\pi b[1exp(N)]𝑑b$$
(2)
where the transport function T(b) can be obtained from the average n-n collision number N as a function of the impact parameter b. More details could be found in Ref. .

In the previous IDQMD model, the volume occupied by nucleon in the projectile and target (Pauli volume) was a constant ($`\mathrm{}^3/2`$). This constant could be deduced from the lowest limit of the uncertain relationship between the momentum and the coordinate. Recently it was found that the Pauli volume is sensitive to the incident energy and should not be a constant qualitatively. Calculation of the average collision number $`N`$ of the $`{}_{}{}^{12}C`$ +$`{}_{}{}^{12}C`$ system at high incident energy shows that the average collision number still has an uptrend after the evolution time of 200fm/c. However, other studies indicated that the average collision number has been saturated after 50fm/c at high incident energy. This indicates that the invariable Pauli volume used in the previous IDQMD model may be not suitable and the Pauli volume may be energy-dependent. From formula (2) we know that $`\sigma _R`$ can be decided by the average collision number $`N`$. It is also clear that N is sensitive to the probability of the Pauli blocking in the collisions. Thus, the energy dependence of the Pauli volume becomes very important in the calculation of $`\sigma _R`$. Up to now, there has not yet been parameterized formula to describe the energy-dependent Pauli volume.

We calculate the cross section $`\sigma _R`$ of the $`{}_{}{}^{12}C`$ \+ $`{}_{}{}^{12}C`$ reaction system in an IDQMD framework in a wide energy range as shown in Fig. 1(a). Here soft EOS and in-medium $`\sigma _{nn}`$ are used. The initialized density distribution of $`{}_{}{}^{12}C`$ comes from the SHF calculation. The solid line with open circles shows the calculated results with the constant Pauli volume. It is obvious that they can not reproduce the experimental data. With adjusted Pauli volume values in the IDQMD calculation we fit the experimental data at different energies \[see the dashed line with open squares in Fig.1(a)\]. From these fits, the Pauli volume values were obtained in a wide energy range as shown in Fig. 1(b). Obviously, the Pauli radius R increases generally at intermediate energy and becomes saturated above 400 MeV/u. By using this energy-dependent Pauli volume, it is found that the average collision number N will be a constant after 50fm/c at high energy and after 100fm/c at intermediate energy, which is consistent with the other studies.

In order to test the modified IDQMD model, we use it to study $`\sigma _R`$ for other reaction systems. In Fig.2, the solid line shows the experimental matter density distribution with a long tail and the dashed line shows the one calculated with the SHF density distribution. The experimental density distribution of $`{}_{}{}^{11}Li`$ is introduced into the IDQMD model to replace the SHF one. Figure 3(a) shows the calculated results of $`\sigma _R`$for the $`{}_{}{}^{11}Li`$ \+ $`{}_{}{}^{12}C`$ reaction system with different density distributions as shown in Fig.2. Clearly, the calculation with the experimental density distribution of $`{}_{}{}^{11}Li`$ can fit the experimental $`\sigma _R`$ better. For the reaction induced by halo nucleus$`{}_{}{}^{11}Li`$, the results with the SHF density distribution are lower than the experimental ones at high incident energy for about a few percent since the SHF calculation of $`{}_{}{}^{11}Li`$ does not give the long tail as experimental one, which is expected to play an important role in the calculation of $`\sigma _R`$. This indicates that the density distribution calculated with the SHF density distribution is not appropriate to those halo nuclei  in the IDQMD calculation. It is obvious that the calculated $`\sigma _R`$ at several hundreds of MeV/u is more sensitive to the long-tail density distribution. Since the central densities are also different between them, more conclusions about the density effect on $`\sigma _R`$ can be obtained by further studies. The research along this line is in progress.

Figure 3(b) shows the calculated results of excitation function of $`\sigma _R`$for the $`{}_{}{}^{12}C`$ \+ $`{}_{}{}^{27}Al`$ reaction system. In Fig. 3, all the filled circles represent the experimental data. The solid line with open squares shows the calculated results with the IDQMD, where the initialized density distribution of $`{}_{}{}^{27}Al`$ and $`{}_{}{}^{12}C`$ are calculated by using the SHF model. Figure 3 shows that the IDQMD calculation with the energy-dependent Pauli volume and the density distribution of the SHF model gives a good trend of excitation function of $`\sigma _R`$. For reaction induced by stable nuclei, the calculated results can fit the experimental values well in a wide energy range. With the experimental density distribution the IDQMD model can also give the fine results of $`\sigma _R`$ induced by halo nuclei.

Figure 4 gives the calculation of $`\sigma _R`$of Li isotopes on the $`{}_{}{}^{12}C`$ target at 790MeV/u. The density distribution of $`{}_{}{}^{11}Li`$ is from the experiment and the others are calculated with the SHF model. It shows that $`\sigma _R`$ varies smoothly with mass number up to 9Li. There is a sudden increase of $`\sigma _R`$ between 9Li and $`{}_{}{}^{11}Li`$ which is corresponding to the halo structure in $`{}_{}{}^{11}Li`$. Our calculation results reproduce both the experimental values and the sudden change between 9Li and $`{}_{}{}^{11}Li`$ fairly well.

In conclusion, the IDQMD model has been introduced to study the total reaction cross section $`\sigma _R`$ by using the energye-dependent Pauli volume which is deduced from the excitation function of $`\sigma _R`$ of $`{}_{}{}^{12}C`$+$`{}_{}{}^{12}C`$. The calculated results can reproduce the experimental $`\sigma _R`$ quite well. It is interesting to investigate further the energy and isospin effects of Pauli volume in collision. For halo nuclei, the calculated results by using experimental density distribution are better than that by using the SHF density distribution. It is suggested that the long tail of the density distribution plays an important role for the halo nuclei and $`\sigma _R`$ at high energy is sensitive to the long tail distribution of halo nuclei. Since the IDQMD model incorporates the isospin-dependences of mean-field, nucleon-nucleon cross section, and Pauli blocking, it is also interesting to study each isospin-effect on $`\sigma _R`$ and other physical quantities. In the IDQMD calculation, the clusters can be judged by the relative momenta and coordinates with an isospin-dependent modified coalescence model and the momentum distribution of the projectile fragments can be obtained easily. Thus, the IDQMD model can be used to study the total reaction cross section, fragmentation cross sections and momentum distribution of fragment of halo nuclei simultaneously. It can be used to give a more comprehensive and reliable criterion of halo structure, which is of great significance. These works are in progress.

One of the authors (Y. B. Wei) gratefully acknowledges Dr. F. S. Zhang for his enthusiastic help in personal discussion.