text:
# Entanglement Simulations of Shor’s Algorithm

## Introduction

Quantum entanglement  is a basic resource in quantum information processing. While its role in quantum communication tasks is quite well understood the same cannot be said about quantum computation. While it is generally believed that entanglement is necessary to achieve an exponential speedup of a quantum algorithm over a classical algorithm  the exact mechanism by which this may happen is unclear. In fact, so far it has not been proven strictly whether entanglement is really necessary for an exponential speedup.

Generally, the argument for the power of quantum computation relies on the assumption that any algorithm that simulates the time evolution of a quantum system will be exponential in the number of quantum bits involved, if it explores the whole state space. The reason is that the dimension of the total state space of a quantum system grows exponentially with the number of subsystems. For pure states, this argument implies that the quantum system needs to evolve into an entangled state to be difficult to be simulated. If it is always in a product state, then the number of parameters required to describe it grows only polynomially with the number of subsystems. For mixed states, however, the situation changes significantly. The set of disentangled states, i.e. the separable states, has the same dimension as the set of all states, although its relative size with respect to the total state space decreases rapidly . Therefore, one could imagine that there are dynamics that always leave the system in a separable state, but which are nevertheless difficult to simulate on a classical computer simply by the fact that the number of parameters that is required to describe the quantum system grows exponentially. As a consequence it may conceivably be possible to have efficient quantum computations on separable states as the algorithm is able to efficiently simulate itself. Furthermore this points to the possibility that efficient quantum computation is possible on mixed states.

Recently Knill and Laflamme  investigated the power of quantum computations using one pure qubit and a supply of maximally mixed qubits and were able to construct a problem that could be solved more efficiently using these resources than any known classical algorithm. Also Schulman and Vazirani  were able to show that given a supply of thermal states one could produce a single pure qubit together with many maximally mixed qubits. The latter could then be discarded and the pure qubit combined with other pure qubits in a quantum algorithm. Indeed NMR quantum computation  does start with initially thermally mixed qubits, although the computation efficiency falls off exponentially with the number of qubits. However, we still have the possibility that these mixed qubits could be used in a useful computation like that of Knill and Laflamme. It is therefore of interest to explore this idea further and see what degree of mixing a quantum computer can tolerate before it loses it’s efficiency.

This paper is an exploration of this idea. We study the efficiency of Shor’s algorithm when the quantum computer is in a highly mixed state. We arrive at the conclusion that Shor’s algorithm can be run on extremely mixed states without significant loss of computational efficiency. Nevertheless, it turns out that despite the significant degree of mixedness Shor’s algorithm runs through some weakly entangled states, leaving the question open as to whether a quantum computer really requires entanglement to be efficient.

The sections of this paper are organized as follows: in section I we give an outline of Shor’s algorithm together with possible gate layouts and interpretations; in section II we examine what is known about simulating quantum algorithms; section III looks at entanglement measures for mixed states and section IV multipartite entanglement and it’s quantification; sections V and VI give details of the simulations used in this work and some of the results obtained and finally section VII gives concluding remarks.

## I Outline of Shor’s Algorithm

### A The Algorithm

Shor’s algorithm for factoring an integer $`N=pq`$, where $`p`$ and $`q`$ are prime, relies on finding the period, $`r`$, of the function $`f_a(x)=a^x\text{mod}N`$, where $`a`$ is some integer less than and coprime to $`N`$ chosen at random. Then, with a sufficiently high probability, at least one of the unknown factors of $`N`$ is given by $`\text{gcd}(a^{r/2}\pm 1,N)`$ which can be calculated efficiently using Euclid’s algorithm. Classically all known algorithms are unable to solve the period finding problem in time polynomial in $`\mathrm{log}N`$, the length of the number being factorized.

The quantum period finding algorithm works as follows: two quantum registers are required whose state spaces are of size at least $`N^2`$ and $`N`$ respectively. As we will be using qubits we will require $`L=2\mathrm{log}_2N`$ and $`n=\mathrm{log}_2N`$ of these two level systems for the two registers respectively. We initially prepare the first register in an equal superposition of all possible states by preparing each of the qubits of the register in state $`\frac{1}{\sqrt{2}}\left(|0+|1\right)`$. The second register is prepared in state $`|1`$, where all the individual qubits are in state $`|0`$ except the first which is in state $`|1`$. We now unitarily transform the two registers with the transformation $`U|x,b|x,ba^x\text{mod}N`$ so that

$$U\frac{1}{\sqrt{t}}\underset{x=0}{\overset{t1}{}}|x,1\frac{1}{\sqrt{t}}\underset{x=0}{\overset{t1}{}}|x,a^x\text{mod}N$$
(1)
where $`t=2^L`$. Now, an inverse quantum Fourier transform

$$F^1|y\frac{1}{\sqrt{t}}\underset{z=0}{\overset{t1}{}}e^{2\pi iyz/t}|z$$
(2)
on the first register of Eq. 1 yields the state

$$\frac{1}{t}\underset{x=0}{\overset{t1}{}}\underset{z=0}{\overset{t1}{}}e^{2\pi ixz/t}|z,a^x\text{mod}N$$
(3)
and the first register now contains information about the period of the function $`f_a(x)`$. We access this information by simply measuring the first register in the number, or computational basis obtaining the result $`|c`$, say. It was then shown in  that the fraction $`c/t`$ is, with a sufficiently high probability, most closely approximated (using the continued fractions method ) by a fraction $`j/k`$ (with $`k<N`$) which in lowest terms has $`k=r`$, the period we are trying to find and will therefore sufficiently often give us a factor of $`N`$.

### B Decomposition into basic gates

These are the essential details of Shor’s algorithm as it was first formulated. We must now, of course, be sure that the algorithm runs in time polynomial in $`\mathrm{log}N`$ for general $`N`$ as well as using polynomial space as has been shown above. The time taken to perform the algorithm is generally assessed by counting the number of basic operations involved. To do this we must have a decomposition of the $`U`$ and $`F^1`$ transformations into elementary gates acting on a small numbers of qubits, usually one, two or three.

One condition on these elementary gates is that each of them is reversible, that is, given the output states we could work out the input states (and obviously vice versa). The condition is imposed by the unitarity, and therefore the reversibility of the transformations $`U`$ and $`F^1`$. Examples of such reversible gates are CNOT’s (2 qubits), TOFFOLI’s (3 qubits) and an infinite variety of 1-qubit gates and 2-qubit gates where a single qubit transformation is ’controlled’ by a second qubit.

Detailed polynomially efficient gate layouts for the modulo exponentiation transformations can be found in . They highlight one other complication: the efficient decomposition of general unitary transformations into small basic gates seems to require auxiliary qubits. These are used during the transformation to store quantum information temporarily but are left in their initial states at the end. Some of these must be prepared in a known state ($`|0`$, say), others may be prepared in a completely unknown, or maximally mixed state which may or may not be entangled to other systems outside the computer. Either way they must be returned to their initial state after use during the computation. If they are not returned to their initial states (or some other known state that is disentangled from the rest of the computer) these qubits will be holding information about the states of the non-auxiliary qubits so that the transformation as viewed on the non-auxiliary qubits alone cannot be unitary or reversible and quantum information will have leaked out of the quantum computer.

A polynomially efficient gate decomposition of the (inverse) Fourier transform into 1-qubit Hadamard transformations and 2-qubit controlled phase rotations can be found in Fig. 1 and is all the gates to the right of, and including, the first Hadamard transform (H). It requires $`O\left((\mathrm{log}N)^2\right)`$ 1- and 2-qubit gates to perform the transformation and is therefore again polynomially efficient in time.

The careful (or experienced) reader will notice that for increasing $`N`$ the conditioned phase rotations are of increasingly small values (controlled $`R_L=\left(\begin{array}{cc}1& 0\\ 0& e^{2\pi i/2^L}\end{array}\right)`$ transformations are used) which would require that the accuracy of the gate implementation is exponential in $`\mathrm{log}N`$. This would require exponential resources (in terms of time/energy etc.) but it clear that the Fourier transform implemented without performing the controlled phase shifts to such high accuracy does not affect the transformation too much and so does not reduce the efficiency too much. In return, however, in a practical situation involving decoherence it is an advantage not to carry out these small phase shifts as the computation will then suffer less errors due to it’s shorter computation time .

### C The phase kickback interpretation

In Fig. 1 the controlled modulo exponentiation of the classical number $`a`$ has been decomposed into $`L`$ successive controlled modulo multiplications by $`a^{2^{L1}}\text{mod}N,a^{2^{L2}}\text{mod}N,\mathrm{},a^2\text{mod}N,a\text{mod}N`$ . We will write these as $`{}_{c}{}^{}U_{a}^{2^{L1}},{}_{c}{}^{}U_{a}^{2^{L2}},\mathrm{},{}_{c}{}^{}U_{a}^{2^1},{}_{c}{}^{}U_{a}^{}`$, where

$`{}_{c}{}^{}U_{a}^{2^x}|0,b`$ $`=`$ $`|0,b`$ (4)
$`{}_{c}{}^{}U_{a}^{2^x}|1,b`$ $`=`$ $`|1,a^{2^x}b\text{mod}N.`$ (5)
The modulo multiplications can be written in this way, as powers of the gate $`U_a`$, because multiplication by $`a^{2^x}`$ modulo $`N`$ is equivalent to multiplying by $`a`$ modulo $`N`$, $`2^x`$ times. Actually performing the modulo multiplications in this way would of course require exponentially many repetitions of the basic gate $`{}_{c}{}^{}U_{a}^{}`$ and is therefore a highly inefficient method but each of the controlled modulo multiplications can be performed in time polynomial in $`\mathrm{log}N`$ after classical precalculation of the numbers $`a^{2^x}\text{mod}N`$ .

Viewing the algorithm as such leads us to another interpretation of the period finding algorithm that is just a change of basis away from Shor’s original formulation, as outlined above. The operation $`U_a`$ has a set of $`r`$ eigenstates $`|\psi _j`$ ($`j=0,\mathrm{},r1`$) with eigenvalues $`e^{2\pi ij/r}`$. Applying a controlled-$`U_a`$ gate to the state

$$\left(|0+|1\right)|\psi _j$$
(6)
(aside form normalization) kicks the acquired phase onto the control qubit:

$$U_a\left(|0+|1\right)|\psi _j=\left(|0+e^{2\pi ij/r}|1\right)|\psi _j.$$
(7)
We cannot prepare the eigenstate $`|\psi _j`$ as this would require knowledge of $`r`$. However, using the result $`_{k=0}^{r1}|\psi _j=|1`$  gives

$${}_{c}{}^{}U_{a}^{}\left(|0+|1\right)|1=\underset{k=0}{\overset{r1}{}}\left(|0+e^{2\pi ij/r}|1\right)|\psi _j.$$
(8)
A measurement of the control qubit, in some chosen basis, will now yield information about the fraction $`j/r`$ for some $`j`$ selected at random, although only one bit of information will be acquired. More information can be acquired about the phase if we perform the controlled- $`U_a^{2^{L1}},U_a^{2^{L2}},\mathrm{},U_a^{2^1},U_a^{2^0}`$ gates, using different control qubits for each, the inverse Fourier transform on the control qubits  and a projective measurement on each control qubit. This will sufficiently often allow us to obtain $`r`$ as described above: by finding the fraction $`j/r`$ closest to $`c/t`$ where c is the result of the measurements on the control qubits.

Shor’s algorithm, then, can be seen in terms of the production and measurement of relative phase information which is related to $`r`$.

### D Using mixed states

In  it was shown that $`U_a`$ also has other sets of eigenstates $`|\psi _{j_d}^d`$ ($`j_d=0,\mathrm{},r_d1`$) with eigenvalues $`e^{2\pi ij_d/r_d}`$ and that in fact nearly all of these (at least $`(p1)(q1)`$ of them) have $`r_d=r`$. Consequently the lower qubits (the 2nd register) in Fig. 1 need not be prepared in the initial state $`_{k=0}^{r1}|\psi _j=|1`$ but can be prepared in the completely unknown state

$$\frac{\mathrm{𝟏}}{N}=\frac{1}{N}\underset{j_d,d}{}|\psi _{j_d}^d\psi _{j_d}^d|$$
(9)
(here we are equating $`|\psi _{j_1}^1=|\psi _j`$). This mixed state algorithm is run exactly as before and the period is found at least $`\frac{(p1)(q1)}{N}`$ times as efficiently as the original pure state algorithm, this factor approaching unity as $`p,q\mathrm{}`$.

This also, in fact, means that any randomly selected state, whether pure or mixed, entangled or not, may be used as an input state for the lower qubits and, on average, the algorithm will run efficiently. These lower qubits may also be mixed because they are entangled to systems outside the computer .

In terms of the number of pure qubits that are needed in the algorithm we should once again address the matter of the decomposition of the controlled-$`U_a`$ transformations into basic gates (Section I B). It is well known that polynomially efficient decompositions of these transformations do exist  but they require auxiliary qubits which it seems need to be prepared in a pure state. With some alterations to these decompositions, however, it has recently been shown that, of those auxiliary qubits that cannot be removed from the algorithm, only one need be prepared in a pure state  (this is not to be confused with the ‘one’ pure qubit of the Abstract and the next section). The rest can be prepared in maximally mixed states, although most can no longer be considered as being auxiliary to the computation.

For simplicity, however, we will not consider these auxiliary qubits or this altered form of the controlled-$`U_a`$ transformations for the rest of this work.

### E Using only one control bit

One further modification to the algorithm is to use a semi-classical Fourier transform  \- notice that the gates of the Fourier transform, both one and two qubit, occur sequentially on the qubits. We could thus replace all of the first register of qubits with a single control qubit and perform the gate operations as follows (see Fig. 2): we implement the first controlled-$`U_a^x`$ gate and the Hadamard transformation and measure the control qubit; after resetting the state of this qubit to $`|0+|1`$ we implement the next controlled-$`U_a^x`$ gate and replace the 2-qubit controlled phase shift with a single qubit phase shift if the result of the first measurement was $`|1`$; we continue in this manner with a Hadamard transformation, single qubit phase shifts given the results of all previous measurements, and another measurement and resetting. At the end of the algorithm we have a set of measurement results that have an identical probability distribution to the algorithm using $`L`$ control qubits in the first register, as in Fig. 1.

## II Simulating algorithms

We saw in the previous section how a quantum algorithm consisting of quantum gates on quantum systems can factorize an integer $`N`$ in time and number of qubits polynomial in $`\mathrm{log}N`$. The question then arises as to whether or not we can turn this into a classical algorithm by writing out the effect of the gates and measurements on a classical system such as a computer or piece of paper. If we then find that we can do this efficiently (in time polynomial in the number of qubits) then we have an efficient classical algorithm for factoring. As no efficient algorithm is known we might be fairly certain that no such efficient simulation is possible and that all simulations of Shor’s algorithm will be difficult.

Of course we already have a ’simulation’ of Shor’s algorithm as we have written down the equations (1) and (3) but we cannot in general derive any of the properties of the computer, or indeed the probability distribution of the measurement results, without writing out the density matrix for the whole system (or doing something else equivalently difficult) at the relevant point.

We would like to look at the system after certain gates or sets of gates. We can easily write out the effect of a single-qubit gate $`A`$ on a single qubit, whether pure or mixed, by writing down the unitary matrix for $`A`$ and working out its effect on the state vector $`𝐯`$ or density matrix $`\rho `$:

$`𝐯^{}`$ $`=`$ $`A𝐯`$ (10)
$`\rho ^{}`$ $`=`$ $`A\rho A^{}.`$ (11)
However, we can only do this if the qubit is completely disentangled from other systems which we may later wish to use for information processing. This is because we clearly cannot, for example, simulate the operation of two single qubit gates on two qubits in an entangled state by tracing out the opposite qubit, implementing the gates on each and taking the combined state after the operations as the tensor product of the two states (this is not even the case if the two operations are the identity). That is, there are entries in the density matrix of the combined system which refer to the combined system rather than to the individual systems themselves.

Consider some examples in Shor’s algorithm. The single-qubit gates in the Fourier transform are likely to be acting on qubits which are entangled to other qubits in the computer (possibly many of them) so to simulate the algorithm correctly we must not consider just the state of the single qubit (which will be mixed in general as it is entangled to other qubits). In contrast we noted in section I D that we may use mixed states in the lower register of Shor’s algorithm and that these may be mixed because they are entangled to systems outside the computer. In this case we can address only those qubits which are within the computer, even though they may be entangled to outside systems, because we will not later be concerned with these outside systems.

Let us consider the problem in more detail, considering pure states first. The state vector of an $`M+1`$-qubit system in a general entangled state requires $`O\left(2^{M+1}\right)`$ complex numbers to be written down and stored . If we wish to classically simulate the application of the single-qubit gate $`A`$ to a qubit that is entangled to $`M`$ other qubits one method of simulation would be to apply the 2 by 2 matrix representing the gate to each of the $`M`$ pairs of amplitudes in state vector corresponding to different states of the remaining $`M`$ qubits. This involves $`M`$ 2 by 2 matrix multiplications and therefore requires $`O\left(2^M\right)`$ operations

So to just simulate the effect of a single-qubit transformation in general takes classical resources exponential in $`\mathrm{log}N`$, if entanglement exists across the $`O\left(\mathrm{log}N\right)`$ qubits.

Of course this argument also follows for mixed states - we cannot simulate the effect of even one single-qubit gate efficiently if it is entangled to many other qubits within the computer. Using a method similar to that described above we can perform the 2 by 2 matrix multiplication on $`\left(2^M\right)^2`$ blocks (pre-multiplication by $`A`$ and post-multiplication by $`A^{}`$) of the density matrix for the whole computer, thereby requiring $`O\left(2^{2M}\right)`$ operations.

But the situation here is more complicated: for pure states it is relatively easy to see if entanglement exists in a simulated algorithm (although it is by no means trivial) and whether therefore the above general method needs to be used to simulate gates acting just on parts of the computer. For mixed states, however, it is harder to determine which qubits are separable. For two qubits to be separable there must exist a decomposition into pure states where the pure states are all separable:

$`\rho _{12}\text{separable}`$ $``$ $`|\psi _i_{12},p_i(>0)\text{such that}`$ (12)
$`\rho _{12}`$ $`=`$ $`{\displaystyle \underset{i}{}}p_i\left(|\psi _i\psi _i|\right)_{12}\text{ where}`$ (13)
$`i|\psi _i_{12}`$ $`=`$ $`|\psi _i^1_1|\psi _i^2_2`$ (14)
For almost all states there will be a decomposition into non-separable states. Finding a separable decomposition however, and indeed finding if it exists, is a difficult task.

A mixed state algorithm, then, may be seen as a mixture of pure state algorithms and even if at each stage these pure state algorithms are entangled there may exist other sets of pure state algorithms which are not entangled which mix to give the same mixed state algorithm. And these disentangled sets of pure algorithms will be different after each gate so we cannot easily preclude one existing .

## III Entanglement measures and mixed states

### A Entanglement measures for pure states

Let us first deal with how we would measure entanglement between two quantum systems whose combined state is pure. Many entanglement measures in some way view entanglement as a resource. The process of teleportation  is an important example of a phenomenon that requires entanglement to be observed at all, as it is required in many other applications of quantum information processing including entanglement swapping, dense coding, precision measurements and hiding classical information  and in the violation of Bell’s inequalities .

To do perfect teleportation of an unknown single qubit state requires one of the four Bell states (or ’EPR pairs’) (Eq. 15 and 16) to be shared between the two separated parties 1 and 2:

$`|\varphi _{12}^\pm ={\displaystyle \frac{1}{\sqrt{2}}}\left(|0_1|0_2\pm |1_1|1_2\right)`$ (15)
$`|\psi _{12}^\pm ={\displaystyle \frac{1}{\sqrt{2}}}\left(|0_1|1_2\pm |1_1|0_2\right).`$ (16)
Any other state of two qubits (which cannot be transformed into one of the above by unitary transformations performed by the two parties separately) cannot be used to perform teleportation perfectly.

To quantify the entanglement of a general state $`|\psi _{12}`$ we could look at how well they perform the teleportation process (with some sort of ’fidelity’ measure between the input state and the output state or the maximum probability for perfect teleportation). In fact, we would rather ask how many Bell states can be obtained from the given state $`|\psi _{12}`$ using only Local quantum Operations (such as transformations, addition and removal of local separable systems and measurements) and Classical Communication (LOCC for short) . The Bell states, therefore, have entanglement of value ’1’ because you can only obtain one Bell state from each Bell state supplied (you cannot on average increase the number of Bell states with LOCC).

For other states of two qubits, say

$$|\psi _{12}=a|00_{12}+b|11_{12},a,b^+,a>b,$$
(17)
(with $`a^2+b^2=1`$) let us consider first what we can do with just one copy of the state. The most efficient method of obtaining a Bell state from this state is to use the Procrustean method . This only creates a Bell state with probability $`2|b|^2`$, the rest of the time creating a completely separable state.

We can do better than this efficiency if we are supplied with more copies, $`n`$ say, of the state $`|\psi _{12}`$ held by the separated parties. Now we can increase the number of Bell states obtained between the parties per copy of $`|\psi _{12}`$ held by allowing each of the separated parties to perform joint operations on those parts of the entangled states each holds. So say $`k`$ is the number of Bell pairs obtained, then on average, and with a suitable method, $`k/n2b^2`$. The exact results for any finite number of copies are known  and asymptotically, in the limit of large n (provided we allow an arbitrarily small probability of error), the average number of Bell pairs obtained per copy for pure states is given by

$$\frac{k}{n}=E\left(\rho _{12}\right)=S\left(\rho _1\right)=S\left(\rho _2\right)$$
(18)
where $`\rho _{12}=\left(|\psi \psi |\right)_{12}`$, $`S(\sigma )`$ is the von Neumann entropy of the density matrix $`\sigma `$, and $`\rho _1`$ and $`\rho _2`$ are the partial density matrices of the first and second of the entangled particles. For the state of Eq. (17) this gives $`E\left(\left(|\psi \psi |\right)_{12}\right)=a^2\mathrm{log}a^2b^2\mathrm{log}b^2>2b^2`$. These, and future results also hold for two party systems composed of individual systems of more than two levels.

In the asymptotic limit (only) the process is also true in reverse: given $`k`$ Bell pairs and taking $`k\mathrm{}`$ we can create from them, using LOCC, $`k/E\left(\psi _{12}\right)=n`$ copies of the state $`|\psi _{12}`$. We say, then, that for pure states the asymptotic entanglement of formation, $`E_F`$ (the number of Bell states we need to pay per copy of state $`|\psi _{12}`$ we get in return) , is the same as the asymptotic distillable entanglement $`E_D`$ (the number of Bell pairs we can distill out of the state $`|\psi _{12}`$ per copy of $`|\psi _{12}`$ paid).

### B Mixed states

For mixed states the situation is far less clear. $`E_F`$ and $`E_D`$ can be defined in the same way but for general mixed states they are not equal - you cannot in general get as many Bell pairs out of a state as you would put in (you certainly cannot obtain more or you would be able to locally create entanglement ad infinitum). This is due to the fact that a mixed state is to some extent unknown and the randomness inserted on the formation of the state from Bell states cannot be eliminated unless extra information about the state is obtained.

Another problem with $`E_F`$ and $`E_D`$ defined in this way is that at present there are no analytical methods for calculating either for general mixed states. The only example where an analytical expression does exist for more than some specific subclasses of states is the entanglement of formation of a single two-qubit system .

These entanglement measures are not the only possibilities we could produce. Others exist such as the relative entropy . So what are the conditions for a mathematical object to be called an entanglement measure? One set of conditions that are generally accepted as being sensible for an entanglement measure of a state $`\rho _{12}`$ are as follows :

(i) $`E\left(\rho _{12}\right)=0`$ if $`\rho _{12}`$ is separable.

(ii) Local unitary transformations on $`\rho _{12}`$ leave $`E\left(\rho _{12}\right)`$ invariant, i.e. $`E\left(\rho _{12}\right)=E\left(U_1U_2\rho _{12}U_1^{}U_2^{}\right)`$.

(iii) $`E\left(\rho _{12}\right)`$ cannot, on average, increase under local operations and classical communication.

We may also wish to add one more condition to this list. This says that

(iv) the entanglement measure should be equal to the pure state entanglement measure (Eq. (18)) for all pure states.

Of course not all entanglement measures need obey this condition, indeed the one we will be using below does not. There are many candidates for entanglement measures and all these measures will not agree with each other for mixed states. More importantly, the measures will put a different order on the states , that is, according to one measure of entanglement $`E_1`$ the state $`\rho _{12}`$ may be more entangled that $`\sigma _{12}`$ but according to another measure of entanglement $`E_2`$ the state $`\sigma _{12}`$ could be more entangled than $`\rho _{12}`$:

$`E_1\left(\rho _{12}\right)>E_1\left(\sigma _{12}\right)`$ (19)
$`E_2\left(\sigma _{12}\right)>E_2\left(\rho _{12}\right).`$ (20)
Indeed, it was shown in  that any two entanglement measures that agree for pure states but not for mixed states must put a different order on the states.

We must accept, then, that entanglement measures may put different orders on states. The only alternative is to declare one entanglement measure as the ’correct’ one, which immediately prevents us from examining how we might prepare entanglement and how we might use it.

### C A measure of entanglement for mixed states

For the present work we need a measure that is easy to calculate and has an analytical form for general mixed states. We will call this measure the logarithmic negativity $`E_{neg}`$ . It is defined as follows: first denote the matrix elements of the density matrix $`\rho _{12}`$ in some tensor product basis by

$$\rho _{12}^{ij,kl}={}_{1}{}^{}i\left|{}_{2}{}^{}j\left|\rho _{12}\right|k_{1}^{}\right|l_{2}^{}.$$
(21)
The partial transpose  with respect to system 2 is then defined in this notation as

$$\left(\rho _{12}^{T_2}\right)^{ij,kl}=\rho _{12}^{il,kj}$$
(22)
(the labels $`l`$ and $`j`$ have swapped places). It has been shown  that the positivity of this new matrix (that is, the positivity of all its eigenvalues) is a necessary condition for the state to be separable. Therefore, $`E_{neg}`$ is now defined as the log of the sum of the absolute values of the eigenvalues of the new matrix $`\rho _{12}^{T_2}`$ or $`\rho _{12}^{T_1}`$ (the eigenvalues and therefore this measure are also independent of the particular tensor product basis in which the state is considered). This can be written in more compact form as

$$E_{neg}=\mathrm{log}Tr|\rho _{12}^{T_2}|=\mathrm{log}Tr|\rho _{12}^{T_1}|.$$
(23)
As mentioned above this measure does not agree with the pure state entanglement measure of Eq. (18). However, one particularly useful property of $`E_{neg}`$ is that it is an upper bound for $`E_D`$:

$$E_{neg}E_D.$$
(24)
What is more important is that if $`E_{neg}=0`$ we can be sure that the two party system does not have distillable entanglement, although it is not known whether the reverse statement is true or not: we cannot say that distillable entanglement does exist if $`E_{neg}0`$. Also the fact that $`E_{neg}=0`$ does not mean that the state is separable (there exist states with $`E_{neg}=0`$ that are inseparable, which are known as bound entangled states  as the entanglement cannot be distilled into Bell states but is somehow bound from us).

## IV Multipartite entanglement

Entanglement does not only exist between two-party (bipartite) systems, it can also exist between three or more parties. One example is the three-party GHZ state :

$$|\psi _{GHZ}_{123}=\frac{1}{\sqrt{2}}\left(|0_1|0_2|0_3+|1_1|1_2|1_3\right).$$
(25)
One particular property of this state is that it has no bipartite entanglement in the sense that if we trace out the third party, say, the remaining bipartite state is given by

$`Tr_3\left(|\psi _{GHZ}\psi _{GHZ}|\right)_{123}`$ (26)
$`={\displaystyle \frac{1}{2}}\left(|0_10||0_20|+|1_11||1_21|\right).`$ (27)
The entanglement in this state is of different nature to the one contained in an EPR state because it is impossible to inter-convert GHZ states and EPR pairs reversibly . While one can create a GHZ state from two EPR pairs, one can only ever obtain one EPR pair from one GHZ state. What is not known for three party systems is whether GHZ states and EPR pairs are the only different kinds of entanglement. What is known is that there are in fact more types of multipartite entanglement for systems of more than three parties . Three-party entanglement, then, may contain two-party entanglement and under certain conditions we can locally reversibly (under LOCC) transform between the two .

How then do we measure the entanglement of a multipartite system? Because of the different types of entanglement involved having one general measure for all these types is difficult (the relative entropy suitably redefined for multipartite systems is perhaps one exception  although there are still many problems).

Our approach will be to use a bipartite entanglement measure to verify the existence of entanglement between all the different qubits in the computer. Let us suppose that we have an entanglement measure for bipartite systems, or a way of verifying whether a state is separable or not between the two parts of the bipartite system. We now use this measure on all possible bipartite partitionings of an $`nparty`$ system. For example, for a system of 4 qubits labeled as 1, 2, 3, and 4 the possible bipartite partitionings are

$`1/234,\mathrm{\hspace{0.17em}\hspace{0.17em}2}/134,\mathrm{\hspace{0.17em}\hspace{0.17em}3}/124,\mathrm{\hspace{0.17em}\hspace{0.17em}4}/123,`$ (28)
$`12/34,\mathrm{\hspace{0.17em}\hspace{0.17em}13}/24,\mathrm{\hspace{0.17em}\hspace{0.17em}14}/23`$ (29)
where $`1/234`$ means system 1 is considered as being partitioned from systems 2, 3 and 4. For general $`n`$ there are $`2^{n1}1`$ such partitionings.

Clearly, if there is entanglement between the two sides for at least one partitioning, then the state is entangled. The question remains as to whether the converse is true, that is, if all bipartite partitionings are separable is the state completely separable i.e. the multipartite state can be written as a mixture of product states? For pure states it is clear that this is also true (in fact, you need only look at some particular subset of the possible partitionings). But for mixed states it is not true - there are states that are separable across all bipartite partitionings (and therefore have positive partial transpose (PPT)) but that are not completely separable . However this entanglement cannot be distillable entanglement \- if we could distill it into pure multipartite entanglement by local operations it could in turn be changed into bipartite entanglement between qubits.

Our method in the simulations, then, will be to use the logarithmic negativity measure across all bipartite partitionings of the qubits in the algorithm. This will tell us if any distillable bipartite entanglement exists in the algorithm but will also show us if distillable multipartite entanglement of any form exists. This follows from the fact that this measure in effect verifies whether the state is PPT and is therefore not distillable across any bipartite partitioning. If this is the case then no multipartite distillable entanglement can exist (we cannot distill the entanglement into pure state entanglement) otherwise we would again be able to distill this into (pure) bipartite entanglement of some form.

So, we can indeed verify in this way whether any distillable entanglement of any form exists although we cannot preclude that there exist non-distillable entangled states.

## V The Simulations

### A The Basic simulations

Let us first introduce the basic method of our simulations. The states and gates are all stored in matrix form, the former as density matrices (because we will in general be using mixed states) and the gates as specific unitary transformations. $`n`$ qubits therefore require $`2^n\times 2^n`$ density matrices and unitary transformations. After inputting the initial state of the computer we simulate the effect of each gate by pre- and post-multiplying the density matrix by the unitary transformation and its Hermitian conjugate representing the gate to obtain the new state of the computer.

#### 1 Simulation of quantum gates

We will not simulate the effect of all single-, two- or three-qubit gates but will only simulate the algorithm as it appears in Figure 2 where we have convenient points at which to examine the computer i.e. we will only simulate the controlled modulo multiplication gates as a single gate and the gates of the Fourier transform. We will also be using this version of the algorithm as it reduces the number of qubits needed (by about 2/3) which will result in a great decrease in time and space resources the algorithm requires to be simulated.

It was noted in section II that, in general, because of the potential entanglement across all the qubits in the algorithm, the simulation of even a single qubit gate requires an exponential number of operations. For example, if the single qubit Hadamard transform

$$H=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)$$
(30)
is acting on the first qubit in a computer consisting of 2 qubits then the unitary transformation required is

$$H\mathrm{𝟏}=\frac{1}{\sqrt{2}}\left(\begin{array}{cccc}1& 1& 0& 0\\ 1& 1& 0& 0\\ 0& 0& 1& 1\\ 0& 0& 1& 1\end{array}\right),$$
(31)
and if it is acting on the second qubit

$$\mathrm{𝟏}H=\frac{1}{\sqrt{2}}\left(\begin{array}{cccc}1& 0& 1& 0\\ 0& 1& 0& 1\\ 1& 0& 1& 0\\ 0& 1& 0& 1\end{array}\right).$$
(32)
To simulate the effect of these gates on the 2 qubit state by straight matrix multiplication however is not optimal. From the form of the above transformations it is clear that it is better to act with the Hadamard transform on blocks. In the case of $`H`$ acting on the first qubit the density matrix written in block form evolves as

$$\rho =\left(\begin{array}{cc}A& B\\ C& D\end{array}\right)\left(\begin{array}{cc}HAH^{}& HBH^{}\\ HCH^{}& HDH^{}\end{array}\right).$$
(33)
Here of course $`C=B^{}`$ (saving us some calculation time) as it is a density matrix and also for the Hadamard transformation $`H=H^{}`$. The same is true for $`H`$ acting on the second qubit except the elements in the 2 by 2 blocks upon which $`H`$ acts (within the 4 by 4 density matrix) are separated as is shown diagrammatically here:

$$\left(\begin{array}{cccc}a& b& a& b\\ c& d& c& d\\ a& b& a& b\\ c& d& c& d\end{array}\right)$$
(34)
(lower case letters correspond to elements of blocks).

For a general $`2^n`$ by $`2^n`$ density matrix the simulation of a single qubit gates thus requires $`O\left(2^{2n}\right)`$ operations. For gates acting on $`mn`$ qubits this number rises to $`O\left(2^{2n+m}\right)`$ operations.

#### 2 Simulation of tracing and measurements

First of all we may need to trace out one (or more) of the qubits (labeled by $`x`$) of the computer to leave the density matrix

$$\rho _{12\mathrm{}(x1)(x+1)\mathrm{}n}=Tr_x\left(\rho _{12\mathrm{}n}\right).$$
(35)
In matrix form this tracing step effectively involves adding two $`2^{n1}`$ by $`2^{n1}`$ sub-matrices together and zeroing the rest of the $`2^n`$ by $`2^n`$ matrix. Thus the tracing step requires $`O\left(\left(2^n\right)^2\right)`$ operations (although these operations are additions rather than multiplications and are therefore considerably quicker).

We can simulate single-qubit projective measurements in a similar way to that of simulating gates. We will assume, without loss of generality, that all measurements are done in the computational basis, as to perform a measurement in a different basis (even an entangled basis) we can unitarily transform the system (with entangling gates if necessary) and measure in the computational basis.

First of all we must calculate the probabilities of a measurement on the qubit yielding the result $`|0`$ ($`p_0`$) or $`|1`$ ($`p_1=1p_0`$). This is done by summing those elements on the diagonal of the density matrix for the whole computer which correspond to the measurement result.

For a Monte Carlo simulation of measurement results we may now generate a random number, $`p`$, from a uniform linear distribution between 0 and 1, and if $`p<p_0`$ take the simulated measurement result to be $`|0`$, otherwise it is $`|1`$.

The measurement changes the state of the computer. To simulate this we must use two projection operators

$$P_0=\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right)\text{and}P_1=\left(\begin{array}{cc}0& 0\\ 0& 1\end{array}\right)$$
(36)
and act on the state with either $`P_0`$ (if the measurement result was $`|0`$) or $`P_1`$ (if the result was $`|1`$) just as we would act with a single qubit gate. This will, of course, result in just setting 3/4 of the element of the density matrix to zero. This gives us the (subnormalised) state of the measurement collapsed computer, including the (partial or total) collapse of any qubits that are entangled to the measured qubit.

We must then be sure to renormalize the state of the whole computer. This can be done easily by dividing each entry of the density matrix for the collapsed computer by the trace of the density matrix (or equivalently the measurement probabilities, $`p_0`$ or $`p_1`$).

The simulation of the measurement therefore takes $`O\left(2^{2n}\right)`$ operations but notice that the quantum computer does this in linear time, as it only needs to find qubit $`x`$ and measure it.

For the full Monte-Carlo type simulation we must of course repeat the simulation of the algorithm a number of times and average any properties of the system we obtain during each simulation.

#### 3 Other processes

We can also re-prepare the measured qubit in the required state using this method by another application of single qubit gates. If we wish to prepare the state $`|0+|1`$ then if the measurement result was $`|0`$ we apply $`H`$ and if it was $`|1`$ we apply a state flip $`\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)`$ followed by $`H`$.

We can now simulate any algorithm we wish with any set of gates and any type of measurement (a Positive Operator Valued Measure (POVM) can be simulated by adding auxiliary systems, performing unitary transformations on these systems together with the computer and measuring the auxiliary systems).

#### 4 Entanglement calculations

We are mainly interested in the degree on entanglement at each stage. As mentioned in section III B the degree of entanglement does not change under local unitary operations so we need only examine the entanglement after operations on more than one qubit. From Fig. 2 this will be after the controlled-$`U_a`$ operations (where entanglement may increase of decrease) as well after the measurements (where the entanglement may also increase or decrease but on average it should never increase - local projective operations can never increase the amount of entanglement on average).

As in section IV we will use the logarithmic negativity measure of section III C, calculated and averaged across all possible bipartite partitionings of the system. If we are using a Monte Carlo simulation these must of course be averaged across repeated simulations of the algorithm as different entanglements will be observed with different measurement results.

#### 5 Efficiency Accountancy

Let us now check the time efficiency (or lack of it) of Shor’s algorithm simulated by the above method:

(i) to form those gates acting on all qubits requires $`O\left(2^{2n}\right)`$ steps and there are $`O(n)`$ of these giving $`O\left(n2^{2n}\right)`$ steps altogether

(ii) for these gates the matrix multiplications for the simulation of the gate require $`O\left(2^{3n}\right)`$ operations. There are $`O\left(n\right)`$ of these gates giving $`O\left(n2^{3n}\right)`$ operations in all.

(iii) each application of a single qubit gate (including measurement projections) takes $`O\left(2^{2n}\right)`$ steps and there are $`O(n)`$ of these, which is $`O\left(n2^{2n}\right)`$ operations altogether.

(iv) at $`O(n)`$ of the steps we wish to examine the entanglement $`O\left(2^n\right)`$ times. This requires $`O\left(2^{2n}\right)`$ steps each for the partial transposition, $`O\left(2^{3n}\right)`$ for the numerical eigenvalue routine  and $`O(n)`$ steps for the remaining calculation of our entanglement measure. From the most inefficient process (numerical eigenvalue calculation) this gives $`O\left(n2^{4n}\right)`$ operations for the whole. This then is the most important stage of the simulation in terms of the efficiency of the algorithm.

### B Tree simulations

In the above accounting we have not included any contribution from the fact that we have to repeat the whole simulation many times for the Monte Carlo method (section V A 2). Let us examine this more carefully. If we wish to estimate the probability distribution $`\{P,1P\}`$ with two possibilities to a certain accuracy, $`ϵ`$, we must repeat the Monte Carlo simulation approximately some number $`T(ϵ)`$ times. At the $`j^{th}`$ measurement of the algorithm there would be $`2^{j1}`$ states the computer could be in at the point (from the previous possible measurement results). For each of these we would need to estimate the probability distribution giving $`2^{j1}\times T(ϵ)`$ simulations we need to perform, this number getting exponentially worse at each step. So, given a fixed number of simulations the accuracy of the probability estimate and therefore the estimate of any properties of the system can get exponentially worse at each measurement step.

This leads us to an alternative simulation method, a ’tree’ type simulation where every possibility of each measurement is considered. We have $`L`$ measurements during the algorithm and at each one the algorithm is given two possible ’paths’ to use for each of the paths already available. So the computer can take two different paths (be in two different states) at the first measurement, four different paths at the second and $`2^L`$ ($`L=2n`$) different paths after the last measurement. This gives us $`2^{2n+1}1`$ possible states of the computer through the algorithm. And we want to sample the entanglement in the computer at two stages between each measurement, after the measurement and after the controlled modulo multiplication gates.

We can calculate the probability of each of the measurement results at each stage of the algorithm and find the probability for each possible path leading to that stage. Having calculated these we can calculate exact averages for any properties of the system at any stage if we know the properties we wish to examine at each branch of the tree.

This gives us, of course, another exponential overhead in our simulation but, as we saw, this was the same for the Monte Carlo simulation.

### C Noise simulation

The high susceptibility of quantum computers to noise from the environment  in which the computer is run has been known for some time. This noise comes from entanglement of the computer with the environment suffered during the execution of the algorithm. Further sources of error are inaccuracies in the measurements and implementation of the gates. Fortunately it has been shown that error-correcting codes exist . These encode the state of a qubit into the joint (entangled) state of many qubits such that random errors occurring independently on the qubits below a certain (reasonable) threshold can be corrected back to the correct state using measurements and unitary transformations based on the measurement results. It has also been shown that these codes can be used in fault tolerant quantum error correcting schemes , that is, the measurements and transformations that implement the error correcting code itself need not be implemented perfectly.

The fact that Shor’s algorithm can be implemented with ’noisy’ mixed states leads us to asking if this in any way increases the algorithm’s robustness to noise, or whether the coherence within the pure state decomposition of the mixed states needs to be preserved.

It will be interesting, then, to simulate the effects of random noise injected into the computer and it’s effect on the efficiency with which the number $`N=pq`$ is factorized.

There are many ways of simulating noise and many types of noise we could use in the simulation. We have selected two types, noise by measurement and noise by random Pauli operations. It should be noted, however, that different types of noise, whether in quantum or classical scenarios, tend to have similar qualitative effects and so we need not worry too much about the precise nature of each.

#### 1 Noise by measurement

Here, we address each qubit in turn and with a given probability apply a measurement on that qubit in the computational basis (although it need not be in this basis for general noise). The probability for the two outcomes is governed by the state of the particle, as described in section V A 2.

This collapses the state of the computer in some way. We will, with the given probability, apply the noise step to every qubit after each gate during the algorithm. Of course, the controlled modulo multiplication gate consists of many gates acting on small numbers of qubits so would have more time to be effected by noise.

We should also carefully note that when running the algorithm we would not know the result of the measurement taken by the measurement noise and the algorithm would at that stage become more mixed (the computer becomes a weighted mixture of the two states post measurement) as we do not know which measurement result was found. In fact we would not even know if a noise step had been applied.

#### 2 Noise by Pauli operations

Here, again addressing each qubit in turn, with a certain probability we will apply one of the three Pauli operators ($`\sigma _x,\sigma _y,\sigma _z`$), or the identity operator ($`I`$):

$`\sigma _x=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)`$ $`\sigma _y=\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right)`$ (41)
$`\sigma _z=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)`$ $`I=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)`$ (46)
to the qubit, these four operations occurring with equal probability.

Again we should be careful to point out that a noise step would leave the computer in a equal mixture of the four states that result after the application of the four Pauli operations. Thus the qubit will have its state completely randomized (although it could still be entangled to other qubits).

### D Mixing of the control qubit

As we want to examine the entanglement in the computer it will be interesting to investigate ways of reducing the entanglement and seeing how this affects the computer. We have already introduced mixed states into part of the quantum computer and seen that it does not affect the efficiency very much. But we could also mix the state of the control qubit. This must affect the efficiency of the algorithm (a totally random algorithm cannot be any use to us) and it will be useful to compare this to the change in entanglement.

We will do this by preparing (and re-preparing) the control qubit in the state

$$(1ϵ)|\psi \psi \left|+ϵ\right|\psi ^{}\psi ^{}|$$
(47)
where $`|\psi =\frac{1}{\sqrt{2}}\left(|0+|1\right)`$ and $`|\psi ^{}=\frac{1}{\sqrt{2}}\left(|0|1\right)`$. Using the methods of our quantum computer simulator this tensor product operation is most conveniently done as follows: we assume that the control qubit is in state $`|0`$ (if it is in the state $`|1`$ after measurement it can be flipped easily); if we denote the state of the whole computer (including the control qubit) by $`\rho _{12\mathrm{}n}`$ we calculate the temporary density matrix

$$\rho _{12\mathrm{}n}^\text{t}=\left(F\mathrm{𝟏}^{n1}\right)\rho _{12\mathrm{}n}\left(F\mathrm{𝟏}^{n1}\right);$$
(48)
where $`F`$ denotes the single qubit flip operator (section V A 3); we can then mix this with the original density matrix for the computer in the required proportions:

$$\rho _{12\mathrm{}n}^{mix}=(1ϵ)\rho _{12\mathrm{}n}+ϵ\rho _{12\mathrm{}n}^\text{t}$$
(49)
and a final Hadamard transformation on the control qubit gives us the required state of the computer.

Changing the parameter $`ϵ`$ between 0 and $`1/2`$ allows us to decrease or increase the mixedness of the computer.

## VI Results

### A Entanglement in Shor’s algorithm

Firstly we will look at the entanglement in the pure and mixed state single control qubit algorithms (calculated using the ’tree’ simulation method (section V B)). We will look at the average bipartite entanglement as measured by the logarithmic negativity entanglement measure (section III C) averaged across all possible bipartite boundaries at each stage of the algorithm (section IV). We will also average across all possible algorithms factoring numbers of 4 binary digits and 5 binary digits. These results are shown in Figs. 5 and 6.

These two sets of results have a very similar form, in particular the entanglement in the mixed state algorithms closely mirrors the mixedness of the quantum computer (the von Neumann entropy of the state of the whole computer) at each stage. Also note that the amount of entanglement increases towards then end of the algorithm where the most significant (i.e. the highest) bits of the number $`c`$ are decided and the least significant bits of the period $`r`$ are found.

We can see immediately that in all algorithms, both pure and mixed, according to our entanglement measure, entanglement exists although it is up to three times lower in Fig. 5 and 6 for the mixed state algorithm. The mixed state algorithms are, however, at least around half as efficient as the pure state algorithms for the values of $`N`$ considered (which can be seen by calculating $`\frac{(p1)(q1)}{pq}`$ for each $`N=pq`$).

### B Noise

Next we examine the effect of noise on particular algorithms, namely for $`N=15`$, $`a=2`$ (Figs. 7 and 8) which has period $`r=4`$ and $`N=21`$, $`a=2`$ (Figs. 9 and 10) which has period $`r=6`$. Both measurement and Pauli noise were implemented as described in section V C. Additionally we considered the situation where the noise was not allowed to act on the control qubit.

For the $`N=15`$ case when noise is allowed to act on all qubits we see an exponential drop off in the efficiency of the algorithm with increasing noise level for both types of noise. However, when noise is not allowed to act on the control qubits the efficiency only falls to the efficiency level of the mixed state algorithm which, as has been noted, is efficient enough. The is due to the way the noise is implemented and the fact that the period is 4. For algorithms with period of the form $`r=2^m`$, for some integer $`m`$, each of the controlled modulo multiplication gates can be applied on any state independently of the state of the system up to that point.

The reason is as follows: when the period is of the form $`2^m`$ the first $`nm`$ controlled modulo multiplication gates are in fact the identity operation (note, however, that individual operations within the gate decomposition of the controlled modulo multiplications will not be the identity operation so noise would then greatly affect the efficiency of the algorithm) and the measurement result after each will be $`|0`$. The next gate is now the only relevant gate. If the measurement result after this gate is $`|1`$ the period will be found correctly. In this case

$$c=\underset{m1\text{ times}}{\underset{}{\mathrm{?}\mathrm{?}\mathrm{?}\mathrm{?}\mathrm{?}\mathrm{?}\mathrm{?}\mathrm{?}\mathrm{?}\mathrm{?}}}1\underset{nm\text{ times}}{\underset{}{000\mathrm{}\mathrm{}000}}.$$
(50)
where ? = 0 or 1. Now, the fraction $`c/t`$ will have denominator $`r`$, independent of the results of these $`m1`$ remaining measurements.

Noise in the pure state algorithm applied to all but the control qubit, then, prepares these qubits in some random pure state (i.e. a mixed state). The first $`nm`$ stages of the algorithm do nothing but we will correctly find the period if the next measurement result is $`|1`$, which will occur with the same probability as for the mixed state algorithm.

This is not the case for algorithms of other periods of course as the initial operations are not identity operations. An example is the $`N=21`$, $`a=2`$ algorithm for which the effects are noise are shown in Figs. 9 and 10

These results show us that for general algorithms we cannot increase the mixing of algorithms during its running by re-preparing another maximally mixed state in the lower register of qubits after each measurement without reducing the efficiency of the algorithm considerably - the previous measurements have prepared a state which must not be significantly altered before the next stage. If the period is of the form $`r=2^m`$ we can do this (although having a period of this form is presumably exponentially unlikely for increasing $`n=\mathrm{log}N`$) but this will have no effect on the entanglement until the last $`m`$ steps. For other algorithms this will reduce the entanglement but will also dramatically reduce the efficiency.

### C Mixing of the control qubit

Let us now look at the effect of mixing the control qubit, as described in section V D. We will do this for both the ’pure’ state algorithm (where all other qubits are initially pure) and the mixed state algorithm. Figs. 11 and 12 show the average bipartite entanglement throughout the entire algorithm versus the mixing parameter $`ϵ`$ for the pure and mixed state algorithms with $`N=15`$ and $`a=2`$ respectively. Also shown is the probability that each algorithm correctly finds the period, $`r`$.

Notice in particular how the entanglement in the mixed state algorithm approaches zero before the algorithm becomes entirely random. For the pure state algorithm it does not do so until $`ϵ=0.5`$ is reached. The point at which the average entanglement is zero (to machine precision) in the mixed state algorithm is around $`ϵ=0.396`$ for the $`N=15`$, $`a=2`$ case. This illustrates how the randomness in mixed states can completely mask the distillable entanglement even before the algorithm becomes entirely random. For the pure state case we have mixed two pure state algorithms (with orthogonal control qubit states), both of which do produce entanglement, and found that the entanglement does not disappear until we mix them maximally. For the mixed state algorithm we have also mixed two orthogonal algorithms which both contain entanglement but the distillable entanglement has been lost before we lose all information about which algorithm is running.

Compare these with similar results for the algorithms with $`N=21`$ and $`a=2`$ in Figs. 13 and 14.

Again we see that the entanglement in the mixed state algorithm is zero before the control qubit is maximally mixed, although this point occurs at the higher value of around $`ϵ=0.470`$.

Finally we combine the results above to plot the probability of finding $`r`$ against the average entanglement. This is shown in Figs. 15 and 16.

For each particular algorithm, then, we see a definite trend of an increase in entanglement giving an increase in the probability of correctly finding the period, although initially increasing the entanglement from zero does not produce such a large increase in this probability. We also see that the entanglement required for a given probability is lower for the mixed state algorithm. The maximal probability of correctly finding the period is of course lower for the mixed state algorithm but in the limits $`p,q\mathrm{}`$ this is negligibly so.

## VII Conclusions

It has been shown that it is possible to efficiently factorize with Shor’s algorithm using only one initially pure qubit and a supply of initially maximally mixed qubits. We have also seen that for algorithms with small numbers of qubits the mixing of the algorithm remains high throughout. It is then a natural question to see whether any entanglement is involved in the execution of the quantum algorithm. We find that the high degree of mixing, however, does not preclude the existence of entanglement and indeed in these example algorithms entanglement does appear to exist, even when the state of the computer starts and remains in a highly mixed state.

Conversely, if we try to reduce this entanglement by introducing further mixing into the control qubit we do reduce the entanglement in the computer but at the expense of a reduction in efficiency of the computation. The mixing of the control qubit also sheds light on the nature of entanglement itself, that is, we do not need to lose all information about the nature of an entangled state before we are completely unable to extract entanglement from it, as we have seen when mixing two orthogonal entangled mixed state algorithms.

We have also seen that Shor’s algorithm operated on mixed states is nevertheless susceptible to noise of different kinds and the algorithm does not in general appear to have any increased robustness to noise (except where the noise model on particular algorithms is too simplified to be accurate) even though this algorithm can be run using highly mixed states .

What remains an open question is as to how the above effects behave for algorithms of increasing numbers of qubits. Presumably the entanglement does remain for algorithms of large numbers of qubits and the algorithm remains highly susceptible to noise but because of the exponential nature of simulating the algorithms verifying this is an extremely difficult task.

## Acknowledgments

The authors would like to thank R. Jozsa, S. Virmani, V. Kendon and W. J. Munro for helpful comments. This work is supported by the United Kingdom Engineering and Physical Sciences Research Council EPSRC, the Leverhulme Trust, two EU TMR-networks ERB 4061PL95-1412 and ERB FMRXCT96-0066 and the EU project EQUIP.


format: text

text:
# Direct measurement of the ϕ(1020) leptonic branching ratio

## Abstract

The process $`e^+e^{}\mu ^+\mu ^{}`$ has been studied by SND detector at VEPP-2M $`e^+e^{}`$ collider in the $`\varphi (1020)`$-resonance energy region. The measured effective $`\varphi `$ meson leptonic branching ratio: $`B(\varphi l^+l^{})\sqrt{B(\varphi e^+e^{})B(\varphi \mu ^+\mu ^{})}=(2.89\pm 0.10\pm 0.06)10^4`$ agrees well with the PDG value $`B(\varphi e^+e^{})=(2.91\pm 0.07)10^4`$ confirming $`\mu `$$`e`$ universality. Without additional assumption of $`\mu `$$`e`$ universality the branching ratio $`B(\varphi \mu ^+\mu ^{})=(2.87\pm 0.20\pm 0.14)10^4`$ was obtained.

Truly neutral vector mesons play an important role in hadron physics due to their direct coupling to photons. This phenomenon is the basis of the phenomenological Vector Meson Dominance model which successfully describes electromagnetic interactions of hadrons. The key parameters of this model are $`V`$$`\gamma `$ coupling constants. They can be extracted from the vector meson leptonic widths under the assumption that leptonic decay proceeds via one-photon annihilation of the quark-antiquark pair constituting the meson. Leptonic widths also determine the total production cross sections of vector mesons in $`e^+e^{}`$ annihilation and are important for calculation of the hadronic contribution to the photon vacuum polarization .

The $`V`$$`\gamma `$ coupling constant is just one number per vector meson. Could these numbers tell us something non-trivial about the underlying QCD dynamics? Shortly after the 1974 “charm revolution”, Yennie noticed that independently of the vector meson flavor content the following relation holds

$$\mathrm{\Gamma }(Ve^+e^{})/<e_q>^212\mathrm{keV},$$
(1)
where $`<e_q>`$ is the mean electric charge of the valence quarks inside the vector meson $`V`$ in the units of an electron charge. For $`\rho ,\omega `$ and $`\varphi `$ mesons this gives the famous rule: $`\mathrm{\Gamma }(\rho e^+e^{}):\mathrm{\Gamma }(\omega e^+e^{}):\mathrm{\Gamma }(\varphi e^+e^{})=9:1:2,`$ which can be considered as an SU(3) symmetry prediction. The surprising fact here is a relatively high ($`10\%`$) precision of the 9:1:2 rule despite SU(3)-flavor symmetry breaking. Inclusion of charm gives even more badly broken SU(4) symmetry, but Yennie’s relation remains valid with the same precision, which means that SU(4) symmetry still persists for the leptonic widths ratios! Inspired by this strange fact, Gounaris predicted $`\mathrm{\Gamma }(\mathrm{{\rm Y}}e^+e^{})=1.2\text{keV}`$  and was closer to reality than any other author . Current experimental situation with leptonic widths  is shown in Table I.

In the nonrelativistic potential model  the leptonic decay width is given by the Van Royen–Weisskopf formula : $`\mathrm{\Gamma }(Ve^+e^{})=16\pi \alpha ^2<e_q>^2|\mathrm{\Psi }(r=0)|^2/M_V^2.`$

Equation (1) implies then that quarkonium wave function at the origin $`\mathrm{\Psi }(r=0)`$ is proportional to the meson mass $`M_V`$. Note that for Coulomb potential $`|\mathrm{\Psi }(r=0)|^2M_V^3`$, while the linear potential gives $`|\mathrm{\Psi }(r=0)|^2M_V`$. So the leptonic widths tell us that the actual potential appears to be something in between. But even if we postulate such a potential, the relation (1) still has no simple explanation. For light quark systems like $`\rho `$, $`\omega `$, and $`\varphi `$ relativistic corrections are essential. There are also strong interaction corrections governed by the scale dependent $`\alpha _s`$. It was argued  that these corrections modify the Van Royen - Weisskopf formula in the following way:

$$\mathrm{\Gamma }(Ve^+e^{})16\pi \alpha ^2<e_q>^2|\mathrm{\Psi }(r=1/m_q)|^2(10.36\alpha _s(M_V))/M_V^2.$$
(2)
Intuitively, appearance of the constituent quark Compton wavelength $`1/m_q`$ in (2) looks natural, because in relativistic theory a particle cannot be localized within a region smaller than its Compton wavelength . Thus we can expect the quark-antiquark pair to annihilate when approaching each other’s relativistic extents . But this intuitive clarity of (2) doesn’t make an explanation of the remarkable regularity of (1) simpler, because (2) shows that leptonic widths are sensitive to the both nonperturbative and perturbative aspects of QCD. Thus it is not surprising that the leptonic widths become a traditional touchstone for various quark models .

This paper is devoted to the measurement of the leptonic branching ratio of the $`\varphi (1020)`$ meson. There are two leptonic decays: $`\varphi e^+e^{}`$ and $`\varphi \mu ^+\mu ^{}`$. The $`\mu `$$`e`$ universality implies for these decays that $`B(\varphi \mu ^+\mu ^{})=B(\varphi e^+e^{})\times 0.9993`$. Presently only the $`\varphi \mu ^+\mu ^{}`$ decay was measured directly (). There are two PDG values for this decay branching ratio . One of them $`B(\varphi \mu ^+\mu ^{})=(2.5\pm 0.4)10^4`$ is based on the experiments on photoproduction of $`\varphi `$ meson . Another value of the branching ratio $`B(\varphi \mu ^+\mu ^{})=(3.7\pm 0.5)10^4`$ is obtained from $`e^+e^{}`$ experiments . In addition the CMD-2 experiment  has some preliminary results on this decay. One can see that the difference between two PDG values for the decay $`\varphi \mu ^+\mu ^{}`$ is about 2 standard deviations and the accuracy of these results is relatively low. Current branching ratio $`B(\varphi e^+e^{})=(2.91\pm 0.07)10^4`$  is based on measurements of the $`\varphi `$-meson total production cross section in $`e^+e^{}`$ collisions. It was obtained by summation of all $`\varphi `$-meson decay modes: $`\varphi K^+K^{}`$, $`K_SK_L`$, $`3\pi `$, etc.. Up to now the accuracy of $`B(\varphi e^+e^{})`$ was much higher than that of $`B(\varphi \mu ^+\mu ^{})`$. But there is a serious factor limiting the precision of $`B(\varphi e^+e^{})`$ obtained in such an indirect way. It is the interference between $`\varphi `$ meson and other vector states, which description is model dependent. Direct measurement of the $`\varphi e^+e^{}`$ decay in the $`e^+e^{}\varphi e^+e^{}`$ reaction is difficult due to its small probability and huge background from the $`e^+e^{}e^+e^{}`$ Bhabha scattering.

The decay $`\varphi \mu ^+\mu ^{}`$ reveals itself as a wave-like interference pattern in the energy dependence of the $`e^+e^{}\mu ^+\mu ^{}`$ cross section in the region close to the $`\varphi `$-meson peak. The amplitude of the interference wave is proportional to $`B(\varphi l^+l^{})\sqrt{B(\varphi e^+e^{})B(\varphi \mu ^+\mu ^{})}.`$ The accuracy of the $`B(\varphi l^+l^{})`$ measurement in this case is limited only by uncertainty in the calculation of the pure QED part of the $`e^+e^{}\mu ^+\mu ^{}`$ cross section. The 0.2% accuracy claimed in  leads to 0.8% systematic error in the interference amplitude. Large statistics collected by SND detector in the vicinity of the $`\varphi `$ resonance allowed us to make direct measurement of the leptonic branching ratio $`B(\varphi l^+l^{})`$ with the accuracy comparable with that of previous indirect measurements of $`B(\varphi e^+e^{})`$.

Our previous study of the $`e^+e^{}\mu ^+\mu ^{}`$ cross section was done using the 1996 data sample with the total integrated luminosity of $`2.6`$ pb-1 . In 1998 two experimental runs were carried out in the center of mass energy range $`E=9841060`$ MeV in 16 energy points. The collider operated with superconducting wiggler  allowing to increase the average luminosity by a factor of two. Higher luminosity led to relative reduction of the cosmic ray background. The total integrated luminosity $`\mathrm{\Delta }L=8.6`$ pb-1 collected in 1998 corresponds to $`13.210^6`$ produced $`\varphi `$ mesons.

The SND experimental setup is described in detail in ref. . The main part of the SND is a spherical electromagnetic calorimeter. The angles of charged particles are measured by two cylindrical drift chambers (DC). An outer muon system, consisting of streamer tubes and plastic scintillation counters, covers the detector. The integrated luminosity was measured using $`e^+e^{}e^+e^{}`$ events selected in the same acceptance angle as the events of the process under study $`e^+e^{}\mu ^+\mu ^{}`$. The systematic uncertainty of the luminosity measurement is 2%, but its contribution to the systematic error of the interference amplitude estimated using the process $`e^+e^{}\gamma \gamma `$ is only 0.8%.

The primary selection criteria for $`\mu ^+\mu ^{}`$ events were similar to those of our previous work :

* total energy deposition in the calorimeter is more than 270 MeV;
* there are two collinear charged tracks in an event with acollinearity angles in azimuth and polar directions: $`\mathrm{\Delta }\phi <10^{}`$, $`\mathrm{\Delta }\theta <25^{}`$ and with the polar angles within $`45^{}<\theta <135^{}`$ ;
* event is not tagged as $`e^+e^{}e^+e^{}`$ by $`e/\pi `$ separation procedure .

To suppress the background from the processes $`e^+e^{}\pi ^+\pi ^{}`$, $`\pi ^+\pi ^{}\pi ^0`$, $`K_SK_L`$, $`K^+K^{}`$ the outer muon system was used: a requirement for both charged particles to produce hits in muon system renders contribution from this background negligible. For example, the contribution from the process $`e^+e^{}\pi ^+\pi ^{}`$ is about 0.2% in the $`\varphi `$-meson peak.

The cosmic ray background was suppressed by restriction of the time $`\tau `$ measured by outer scintillation counters with respect to the beam collision moment : $`\tau <10`$ ns. About 30% of events selected by the cuts described above are still cosmic ray background. To determine the contribution of cosmic background more accurately the selected events were divided into two classes: (1) $`\mathrm{\Delta }\phi <5^{}`$; (2) $`\mathrm{\Delta }\phi >5^{}`$. The resolution in $`\mathrm{\Delta }\phi `$ is about $`1^{}`$. The $`\mathrm{\Delta }\phi `$ distribution for $`e^+e^{}\mu ^+\mu ^{}`$ events (Fig. 2) was obtained from the experimental data after strong cuts on a difference between time measurements by the muon system for both tracks. Almost all $`\mu ^+\mu ^{}`$ events belong to the first class. The second class contains only $`1.7\%`$ of $`\mu ^+\mu ^{}`$ events. The $`\mathrm{\Delta }\phi `$ distribution for pure cosmic ray events collected in a special run without beams is shown in Fig. 2. The uniformity of this distribution is an artifact of our DC track reconstruction algorithm in which the origin of a charged track in the $`X`$-$`Y`$ plane is fixed to the beam collision point. From Fig. 2 the ratio between numbers of cosmic ray events in the two classes was found $`k_{cs}=N_1^{cs}/N_2^{cs}=1.028\pm 0.033.`$

The number of cosmic ray background events in the first class was calculated for each energy point $`E_i`$ by the following formula: $`N_1^{cs}(E_i)=k_{cs}T(E_i)\mathrm{d}N_2/\mathrm{d}T.`$ Here $`T(E_i)`$ is a data acquisition time for an energy point $`E_i`$, $`\mathrm{d}N_2/\mathrm{d}T`$ is the cosmic event rate in class two averaged over both 1998 experimental runs. The net number of $`\mu ^+\mu ^{}`$ events for each energy point was obtained by subtraction of the cosmic ray background: $`N^\mu (E_i)=N_1(E_i)N_1^{cs}(E_i).`$ The errors of the numbers $`N^\mu (E_i)`$ include the errors of $`N_1(E_i)`$ and $`N_1^{cs}(E_i)`$.

Energy dependence of the detection cross section was fitted according to the following formula:

$`\sigma _{\mu \mu }^{vis}(E)`$ $`=`$ $`\sigma _0(E)R(E)\left|1Z_\mu {\displaystyle \frac{m_\varphi \mathrm{\Gamma }_\varphi }{\mathrm{\Delta }_\varphi (E)}}\right|^2,`$ (3)
$`\sigma _0(E)`$ $`=`$ $`2\pi \alpha ^2\beta (E)(1\beta ^2(E)/3)/E^2,`$ (4)
where $`\alpha `$ is the fine structure constant; $`\beta (E)=(14m_\mu ^2/E^2)^{1/2}`$; $`m_\varphi `$, $`\mathrm{\Gamma }_\varphi `$, $`\mathrm{\Delta }_\varphi (E)=m_\varphi ^2E^2iE\mathrm{\Gamma }(E)`$ are the $`\varphi `$-meson mass, width and inverse propagator respectively; $`\sigma _0(E)`$ is Born cross section of the process $`e^+e^{}\mu ^+\mu ^{}`$; $`Z_\mu Q_\mu e^{i\psi _\mu }`$ — interference parameter. The modulus of the interference parameter is related to the leptonic branching ratio: $`Q_\mu =B(\varphi l^+l^{})3/\alpha .`$ The factor $`R(E)`$ takes into account the detection efficiency and radiative corrections:

$$R(E)=\epsilon _\mu \frac{\sigma _{\mu \mu }(E)}{\sigma _0(E)\left|1Z\frac{m_\varphi \mathrm{\Gamma }_\varphi }{\mathrm{\Delta }_\varphi (E)}\right|^2}.$$
(5)
Here $`\sigma _{\mu \mu }`$ is the result of Monte Carlo integration of the differential cross section of the process $`e^+e^{}\mu ^+\mu ^{}(\gamma )`$  for our geometric cuts with the energy dependent probability for muons to hit outer scintillation counters taken into account. The uncertainty in the energy dependence of this probability adds 1.7% to the systematic error of $`Q_\mu `$. The parameter $`\epsilon _\mu `$ represents energy independent factor in the detection efficiency. It is determined mainly by the cut on total energy deposition in the calorimeter. The value of $`\epsilon _\mu =0.84\pm 0.01`$ was obtained using Monte Carlo simulation of the process $`e^+e^{}\mu ^+\mu ^{}(\gamma )`$ in SND detector , but in the fitting procedure $`\epsilon _\mu `$ was left free. In the calculation of radiative corrections the interference parameter was assumed purely real and equal to $`Z=B(\varphi e^+e^{})3/\alpha =0.120`$ with the PDG table value for $`B(\varphi e^+e^{})`$.

The fitting was performed for two experimental runs independently. Fits with a free $`\psi _\mu `$ yield the interference phase, which is in good agreement with the expected zero value: (1) $`\psi _\mu =(1.0\pm 2.8)^{}`$, (2) $`\psi _\mu =(0.1\pm 2.8)^{}`$. Therefore the interference phase was fixed to $`\psi _\mu =0`$. The fit results presented in Table II show statistical agreement between two experimental runs. Therefore combined fit was performed to obtain the final results which are listed in the third column of the Table II. The values of $`\epsilon _\mu `$ obtained in the fit and from Monte Carlo simulation agree well. The energy dependence of the measured cross section and the fitting curve are shown in Fig. 3. The systematic error of $`Q_\mu `$ includes 1.7% from the uncertainty in the energy dependence of the probability for muons to hit the outer system, 0.8% from the luminosity measurements and 0.8% from the calculation of the radiative corrections. The resulting systematic error is 2%.

In conclusion we obtain the following $`\varphi `$ meson parameters from the measured $`Q_\mu `$ value: $`B(\varphi l^+l^{})=(2.89\pm 0.10\pm 0.06)10^4`$; $`B(\varphi e^+e^{})B(\varphi \mu ^+\mu ^{})=(8.36\pm 0.59\pm 0.37)10^8.`$ This result is in good agreement with our previous one $`B(\varphi l^+l^{})=(3.14\pm 0.22\pm 0.14)10^4`$ . Using PDG value of $`B(\varphi e^+e^{})=(2.91\pm 0.07)10^4`$ we obtain: $`B(\varphi \mu ^+\mu ^{})=(2.87\pm 0.20\pm 0.14)10^4.`$ The good agreement of $`B(\varphi \mu ^+\mu ^{})`$ and $`B(\varphi e^+e^{})`$ confirms the $`\mu `$$`e`$ universality.

This work is supported in part by Russian Foundation for Basic Research, grants No. 99-02-16815 and 99-02-17155. We thank G.V.Fedotovich for useful discussions.


text:
# REFERENCES

ORAL/POSTER REFERENCE: ICF100796OR

FORMAL CONSIDERATIONS ABOUT FRACTURE:

NUCLEATION AND GROWTH

James P. Sethna

Laboratory of Atomic and Solid State Physics,

Clark Hall, Cornell University, Ithaca, NY 14853-2501, USA

ABSTRACT

In this paper, we’ll answer several abstract, formal questions about the nature of crack growth and nucleation. Bringing a field theory point of view to fracture illuminates things in what I hope will be an entertaining way.

Formally, what is the crack nucleation rate? Fracture is an instability of elastic theory under tension. For thermally nucleated cracks, there is an analogy with supercooled liquids. Here the crack nucleation rate can be thought of as an imaginary part of the free energy - giving the decay rate of the metastable, stretched material. As an amusing consequence, we can formally calculate the asymptotic form for the high-order nonlinear elastic coefficients, and explicitly show that elastic theory has zero radius of convergence.

Formally, can we derive the the crack growth laws from symmetry? We can describe mixed-mode three-dimensional fracture as a moving curve in space, decorated with a description of the local crack plane, and driven by the stress intensity factors along the crack. Imposing the symmetries gives us the form of the crack growth laws in two dimensions, explaining the (well-known) fact that cracks under shear (mode II) abruptly turn until the stress becomes purely extensional (mode I). The form of the crack growth law in three dimensions will be derived, and linear stability analysis for mixed-mode fracture will be briefly summarized, with connection to factory-roof morphologies.

KEYWORDS

fracture, crack, nucleation, growth, thermal, evolution, factory roof

INTRODUCTION

In this paper, I will summarize formal work done in collaboration with Jennifer Hodgdon and Alex Buchel on nucleation and growth of cracks in brittle, isotropic, homogeneous materials. This work is formal in the sense that bring to the problem the perspectives and approaches of theoretical physics: we’ll be asking fundamental questions that most of the participants in this conference wouldn’t think to ask. Sometimes answering very basic questions can be illuminating; in any case I hope it will be entertaining.

CRACK NUCLEATION

What is the crack nucleation rate? In a practical context, one may ask what the failure rate is of a component in a machine: the number of breakages per unit time. This definition depends on the size and shape of the object and the external loading. More generally, for a given material under a given load, there will be a rate of crack formation per unit volume per unit time.

In general, the physics underlying this rate is complicated. Cracks always nucleate at inhomogeneities in the material (grain boundaries, surface heterogeneities, pre-existing microcracks, dislocation tangles). These heterogeneities usually evolve in time (fatigue, work hardening, electromigration, void formation). This makes the failure rate history dependent, and too complicated a problem to start with. So, we’ll focus on brittle fracture, where we perhaps can assume that the material far from the crack tip does not irreversibly transform.

With this approximation, we have two cases. (1) The important case is when the failure is due to the inhomogeneities: distributions of microcracks, distributions of bond strengths, or other randomness. In all of these cases, the rate is number of failures for a given increment of strain, not per unit time. Either the stress is large enough to trigger failure or not: a given sample, if it survives the initial loading, will survive forever. This leads to some elegant theories and simulations, but isn’t what I’m focusing on today. (2) The case we study is when the material is homogeneous, and the failure is triggered by a thermal fluctuation. This basically never happens, but it’s an interesting theoretical problem.

So, fundamentally, what is the fracture nucleation rate for our brittle, homogeneous material? This is actually a deep theoretical physics question. The basic problem is defining what is meant at finite temperature by a material under tension. The lowest energy state (free energy, to be pedantic) for a large cube of brittle material under external strain can easily be seen to be broken in two. The elastic energy released is proportional to the volume; the surface energy cost is proportional to the cross-sectional area; for fixed external strain there will be a length at which breaking wins. One should note that this breaking strain in real materials isn’t very large: a cubic meter of iron will prefer to break in two when stretched by 10$`\mu `$m (140 $`\mu `$m if one includes the energy lost due to plastic deformation during fracture).

Crack nucleation is precisely the mechanism by which a material goes from the metastable stretched state to the state which is broken in two. There is a close analogy to the nucleation of crystals in supercooled liquids. In supercooled liquids the liquid is a metastable state with a decay rate given by the thermal probability of generating a critical nucleus. For fracture, the stretched state is a metastable state with a decay rate given by the thermal probability of generating a crack of the size given by the Griffiths threshold.

(The analogy even extends to the practicality of the thermal nucleation mechanisms. Homogeneous nucleation of crystals in supercooled liquids almost never happens: the crystal almost always nucleates on a dust particle, or on the surface…)

Let’s do this calculation explicitly, for a two-dimensional straight crack of length $`\mathrm{}`$, in a medium under external hydrostatic tension (a pressure $`P<0`$). One can solve the elastic equations to find the energy of this crack. Let $`\alpha `$ be the surface energy (energy per unit length, in two dimensions, of the crack surface), $`Y`$ be the Youngs modulus, and $`\sigma `$ be Poisson’s ratio. If we define a critical crack length $`\mathrm{}_c`$ (the Griffiths threshold) by

$$\mathrm{}_c=4Y\alpha /\pi P^2(1\sigma ^2)$$
(1)
then the energy of the crack is

$$E=2\alpha \mathrm{}\alpha \mathrm{}^2/\mathrm{}_c.$$
(2)
It follows that cracks with $`\mathrm{}>\mathrm{}_c`$ will grow to lower their energy, and that cracks with $`\mathrm{}<\mathrm{}_c`$ will heal (at least in our model). The energy at the top of the barrier is $`E(\mathrm{}_c)=\alpha \mathrm{}_c=4Y\alpha ^2/(\pi P^2(1\sigma ^2))`$, and the probability for this critical crack to form by a thermal fluctuation is suppressed by the Boltzmann factor $`\mathrm{exp}(E(\mathrm{}_c)/k_BT)`$. The fracture nucleation rate $`\mathrm{\Gamma }_{failure}`$ per unit volume per unit time will be given by some prefactor times this Boltzmann factor:

$$\mathrm{\Gamma }_{failure}=\mathrm{\Gamma }_0\mathrm{exp}(E(\mathrm{}_c)/k_BT).$$
(3)
(Calculating this prefactor $`\mathrm{\Gamma }_0`$ is quite a subtle problem.)

This calculation is precisely analogous to the corresponding one for a critical nucleus of solid forming in a supercooled liquid. Indeed, the crack as it nucleates will be filled by a gas of vapor atoms (think of fracture in a piece of dry ice, where the vapor pressure isn’t vanishingly small). The spontaneous fracture rate is precisely the nucleation rate of a gas forming in a super-stretched crystal.

Remember, though, that we claimed from a fundamental point of view the difficulty is defining the stretched material. A stretched piece of material is fundamentally a transient: how do we use equilibrium statistical mechanics to study what is only a metastable equilibrium? The precise definition of the stretched, metastable state is as an analytical continuation from the compressed state. In particular, suppose we study the energy $`E(P)`$ of the material under positive hydrostatic pressure $`P`$, where the material is well defined. If we include the effects of thermal fluctuations producing small cracks, and then analytically continue $`E(P)`$ to negative pressure, one can show that it becomes complex! The decay rate (3) is proportional to the imaginary part of the energy:

$$\mathrm{Im}E(P)=C\mathrm{\Gamma }_{failure}(P)=E_0\mathrm{exp}\left(4Y\alpha ^2/\pi P^2(1\sigma ^2)k_BT\right).$$
(4)
Metastable states in statistical mechanics are defined by an analytic continuation of the (free) energy. The nucleation rate is defined as a particular dynamical constant prefactor $`\mathrm{\Gamma }_0/E_0`$ times the imaginary part of this complex energy.

What can we derive from this formal imaginary part to the energy? The real part of the energy is related to the bulk modulus? Let’s define the pressure-dependent bulk modulus $`K(P)`$ by

$$\frac{1}{K(P)}=\frac{1}{V}\left(\frac{V}{P}\right)_T=\frac{1}{PA}\left(\frac{\mathrm{Re}(E)}{P}\right)_T=c_0+c_1P+c_2P^2\mathrm{}+c_nP^n+\mathrm{}$$
(5)
where the elastic material has cross-sectional area $`A`$. $`c_0`$ is the inverse of the normal bulk modulus, and the higher order terms $`c_n`$ represent nonlinear elastic coefficients. These are normally ignored, because by the time they become important plastic deformation or fracture has set in, but in principle a careful experiment could measure them. Now, if $`E(P)`$ is an analytic function, then it obeys Cauchy’s theorem,

$$E(P)=1/2\pi i\underset{\gamma }{}E(T)/T+PdT=1/2\pi i\underset{\gamma }{}E(T)/T\underset{n=0}{\overset{\mathrm{}}{}}(P/T)^ndT.$$
(6)
with $`\gamma `$ circling $`P`$ in the complex plane. Because it divides by $`i`$, it can be used to write the real part of the free energy $`E(P)`$ in terms of the imaginary part. Exchanging the sum and the integral, choosing contours and checking asymptotics, we can calculate the $`n`$-th nonlinear elastic constant in equation 5 asymptotically for large $`n`$:

$$c_n(1)^n(n+2)/\pi A\underset{0}{\overset{B}{}}\mathrm{Im}E(T)/T^{n+1}𝑑T$$
(7)
So because the (complex) energy is an analytic function, the high-order nonlinear elastic constants are related to the fracture rate.

We used this to answer an extremely formal question. What is the radius of convergence of nonlinear elastic theory? We can use the ratio test: the radius of convergence equals the ratio $`c_n/c_{n+1}`$ as $`n\mathrm{}`$. Using (7) and (4) and doing the integrals, we find

$$c_n/c_{n+1}n^{1/2}\left(8\beta Y\alpha ^2/\pi (1\sigma ^2)\right)^{1/2}\text{ as n}\mathrm{}.$$
(8)
(Doing the prefactor carefully with various continuum limits leaves this expression unchanged, except that the surface tension $`\alpha `$ becomes the temperature-dependent surface free energy.) The ratio goes as $`n^{1/2}`$, and the radius of convergence is zero.

Should we be concerned? Hooke’s law is the first term of a series that doesn’t converge? First, the physical origin of this is that our brittle material is in principle unstable as soon as the pressure $`P`$ becomes negative. If you’re not concerned about the fact that your bridges in theory will eventually fall down because of thermal crack nucleation, perhaps the associated non-convergence of Hooke’s law shouldn’t distress you. Indeed, thermal fluctuations are probably too weak to make observable effects in any measurement of the nonlinear elastic constants. Second, nonlinear elastic theory has lots of company: quantum electrodynamics (the most quantitatively successful theory in physics) and Stirling’s series for $`n!`$ also have zero radius of convergence.

Intellectual curiosity aside, what of practical importance can we glean from our calculation? First, perhaps the high-order nonlinear elastic constants of disordered materials might be computed by similar means, and (disorder being important to the fracture rate) might have effects large enough to measure. Second, the imaginary part of the frequency dependent elastic constants are related to viscosities: perhaps considering both external strain and frequency together might be illuminating (vibration-induced thermal fracture rates?) Finally, the corresponding calculation for thermal nucleation of dislocation loops would be straightforward, and might relate to the plastic flow rate at very early stages (before work hardening sets in). The corresponding problem in two dimensions has been well studied, but it seems likely that tangling of dislocations dominates in three dimensions.

CRACK GROWTH LAWS

What are the crack growth laws? On long length and time scales, most macroscopic systems obey simple laws of motion, which depend on the properties of the material only through setting constants (like density and elastic constants) or functions (like surface energies including the crystalline anisotropy). Can we deduce the laws describing the growth of a crack front, coarse-graining over all the atoms, grains, precipitates, and plastic deformation near the crack tip, to get relatively simple laws depending on only a few materials parameters or functions?

Jennifer Hodgdon and I studied this problem, inspired by conversations with Tony Ingraffea and members of his group. Ingraffea’s group solves the hard part of the problem: given a complicated, three-dimensional crack, they can tell us the stress intensity along the crack tip. All they wanted was the rules for advancing the crack forward in time. The simplest problem for us to study, and the one of most current interest to their group, was quasistatic fracture in an isotropic, homogeneous medium. Including crystalline anisotropy would lead us to the physics of cleavage along the crystalline planes. Including random inhomogeneities into theories of brittle fracture leads to logarithmically rough fracture surfaces, but not as rugged as the experiments. Including the inertial dynamics of the growing front leads to interesting traveling waves along the crack front that might improve the agreement with experiment. In ductile and intergranular fracture, the roughness reflects void growth and the presence of grains and inclusions. Our work did not include any of these effects.

Naturally, Ingraffea’s group already had a working solution. For mode I fracture, the material is described by a crack growth velocity $`v(K_I)`$. For mode II fracture, the crack turns abruptly until it is mode I (the principle of local symmetry). There were two other rules in the literature for picking the crack growth direction: one maximizing the energy release, and one moving into the direction of minimum strain energy density. All three rules seemed to predict crack paths which agreed well with one another. At that time (around 1990) they were just moving into 3D fracture, and were quite interested in more crack growth laws for more general geometries.

Jennifer Hodgdon and I decided to try to write the most general crack growth law allowed by symmetry. Let us describe the local crack front with an orthonormal triad of vectors: $`\widehat{𝐛}`$ perpendicular to the local crack face, $`\widehat{𝐭}`$ tangent to the crack front, and $`\widehat{𝐧}`$ pointing in the direction along which the crack was last growing. Our theory is valid on long length and time scales: in particular, we assume that the curvatures of the the radii of curvatures for the crack surface and crack front are large, so to zeroth order the crack is straight and flat (we’ll expand in the curvatures). There are then three symmetries of the crack: the reflection $`R_{\mathrm{plane}}`$ in the $`\widehat{𝐧}\widehat{𝐭}`$ plane of the crack, the reflection $`R_{\mathrm{cross}}`$ in the $`\widehat{𝐧}\widehat{𝐛}`$ plane perpendicular to the crack front $`\widehat{𝐭}`$, and the 180 rotation about $`\widehat{𝐧}`$ (which is the product of the two other symmetries).

The first convenient thing to fall out of the symmetry analysis is the decomposition of the strain field into modes. The linear elastic solution near the crack tip can be decomposed into a superposition of solutions with definite symmetries under $`R_{\mathrm{plane}}`$ and $`R_{\mathrm{cross}}`$. Mode I is symmetric under both transformations, mode II is antisymmetric under reflections in the plane but symmetric to reflections perpendicular to the front, and mode III is antisymmetric under both. Jennifer found that there was yet another mode (mode IV?) which is symmetric under $`R_{\mathrm{plane}}`$ and antisymmetric under $`R_{\mathrm{cross}}`$. Presumably Hodgdon’s mode IV is usually ignored because it does not introduce a crack opening.

We then assume that the body being cracked and the radii of curvature of the local crack front and crack surface are all large compared to the inhomogeneities in the material and the nonlinear zone around the crack tip. This is what is needed to ensure that there is a region around the crack tip where the elastic displacement fields scale with the square root of the distance $`r`$ to the crack tip. (In general the displacement field for each mode can be expanded in a series of all half-integer powers of $`r`$. The less singular terms are important only far from the crack, and are determined by the boundary conditions and the crack shape; the more singular terms are determined by the nonlinear zone around the crack tip, and are also cut off by that zone.) Thus, for example, we describe ductile fracture only if the body is much larger than the plastic zone. We write the growth laws of the crack in terms of the three stress-intensity factors $`K_I`$, $`K_{II}`$, and $`K_{III}`$ in this (so-called $`K`$-dominant) region, their gradients along the crack front, and the curvatures of the crack front and crack surface. We impose the symmetries of the problem described above, plus a gauge symmetry associated with reparameterization of the crack front, to find the most general law of motion allowed by symmetry.

In two dimensions, we found that the most general crack growth law was of the form

$$𝐱/t=v(K_I,K_{II}^2,K_{III}^2)\widehat{𝐧},\widehat{𝐧}/t=f(K_I,K_{II}^2,K_{III}^2)K_{II}\widehat{𝐛},$$
(9)
where $`𝐱`$ is the location of the crack tip, $`v`$ is the velocity of the crack, and $`f`$ represents the tendency for a crack to turn under the external load (hence changing the direction $`\widehat{𝐧}`$ pointing along the crack).

What does this equation imply about the crack growth? If $`K_{II}=0`$, the crack won’t turn; otherwise it will turn (if $`f>0`$) so as to reduce $`K_{II}`$ in magnitude. Thus we agree with the principle of local symmetry, except that our cracks will turn gradually towards this direction, rather than jump abruptly to the new orientation. Over what turning radius does our crack turn? Using results of Cotterell and Rice, Hodgdon showed that if the angle of the crack differs from the angle that makes $`K_{II}=0`$ by a small amount $`\mathrm{\Delta }\theta `$, then $`K_{II}=K_I\mathrm{\Delta }\theta /2`$; this implies that $`\mathrm{\Delta }\theta \mathrm{exp}(fK_Ix/2v)`$, and the angle decays exponentially to pure mode I with a material-dependent decay length of $`2v/fK_I`$. We expect this length scale to be set by a microscopic scale characteristic of the material: the atomic size in a glass, the size of the nonlinear zone in a ductile material, etc. This means, in essence, that our analysis agrees with the traditional prescription that the crack turns abruptly. We think of writing the crack growth laws as differential equations replaces the atomistic short-distance cutoff with a smooth one.

Our analysis did not make any assumptions about microscopic mechanisms. We believe that arguing whether the crack turns to maximize the energy release or minimize the strain energy density is misguided: if you turn the crack according to the wrong rule, in the succeeding step it’ll keep turning and you’ll converge to the same final trajectory (making $`K_{II}=0`$) with a radius of curvature given by the step size in your algorithm.

So much for reinventing what’s already known. What about three dimensions? Now everything depends on where we are along the crack front. If we let $`s`$ measure arclength along the crack front, then up to first order in gradients

$`x/t=`$ $`v\widehat{n}+w\widehat{t}`$ (10)
$`\widehat{n}/t=`$ $`[{\displaystyle \frac{v}{s}}+w{\displaystyle \frac{\widehat{t}}{s}}\widehat{n}]\widehat{t}+[fK_{\mathrm{II}}+g_\mathrm{I}K_{\mathrm{III}}{\displaystyle \frac{K_\mathrm{I}}{s}}+g_{\mathrm{II}}K_{\mathrm{II}}K_{\mathrm{III}}{\displaystyle \frac{K_{\mathrm{II}}}{s}}+g_{\mathrm{III}}{\displaystyle \frac{K_{\mathrm{III}}}{s}}+`$ (11)
$`h_{tb}{\displaystyle \frac{\widehat{t}}{s}}\widehat{b}+h_{nt}K_{\mathrm{II}}{\displaystyle \frac{\widehat{n}}{s}}\widehat{t}+(h_{nb}K_{\mathrm{II}}K_{\mathrm{III}}+w){\displaystyle \frac{\widehat{n}}{s}}\widehat{b}]\widehat{b},`$ (12)
where $`f`$, $`g_\alpha `$, and $`h_{ij}`$ are functions of $`K_\mathrm{I}`$, $`K_{\mathrm{II}}^2`$, and $`K_{\mathrm{III}}^2`$, and the velocity $`v`$ can be a function of these and a number of gradient terms. Only $`f`$ among these constants, which gives the turning radius for mode $`II`$ fracture, is expected to be large (as discussed above): the others should be all of the same order of magnitude.

What can we do with this expression? First, if the crack begins nearly flat and straight, we can do perturbation theory to see if the crack becomes flatter and straighter, or if it goes unstable. Jennifer Hodgdon found that cracks under mode I are stable (as expected), but that cracks under mixed mode I and mode III (twisting) can be stable or unstable depending on the materials constants: in particular, if $`g_I>0`$, steady-state mode III cracks are unstable to small perturbations. Experimentally, she found that cracks under mode III turned to form a helical crack surface: such turning wasn’t allowed for in our perturbation theory, which assumed periodic boundary conditions. Others find, for longer cracks, that mode III fracture is unstable to the formation of a “factory roof” morphology. It wasn’t possible at the time to test Hodgdon’s results with realistic crack-growth simulations, and this work remains unpublished.

Second, we can use it to propagate real 3D cracks in simulations, and compare to experiments and more microscopic simulations. I anticipate that our symmetry analysis will prove important and useful, but that there will be additional internal variables (perhaps crack-tip curvature) that will be necessary for describing the crack, and which may induce crack branching and other effects missing in our existing theory.