text:
# Teleportation via thermally entangled state of a three-qubit Heisenberg XX ring

## Case I: $`B_m=0`$

$$\chi _{ABC}=\frac{1}{Z}\left[|000000|+e^{\beta J}(|W^2W^2|+|W^3W^3|+|W^5W^5|+|W^6W^6|)+|111111|\right]$$

$$\frac{1}{4}(|W^2W^2|+|W^3W^3|+|W^5W^5|+|W^6W^6|)$$
(5)
with $`Z=2+4e^{\beta J}`$. From Eq.(4), the above equally-weighted mixture has

$$C(\rho _{AB})=\mathrm{max}\{\frac{1+e^{2\beta J}}{1+e^{2\beta J}+2e^{\beta J}},0\}=0.$$
(6)
## Case II: $`0<B_m2J`$

$$\chi _{ABC}=\frac{1}{Z}[e^{\frac{1}{2}\beta (B_m2J)}(|W^2W^2|+|W^3W^3|)$$

$$+e^{\frac{1}{2}\beta (B_m+2J)}(|W^5W^5|+|W^6W^6|)+e^{\frac{3}{2}\beta B_m}|111111|]$$

$$=\frac{1}{e^{\beta (2B_mJ)}+2e^{\beta B_m}+2}(|W^2W^2|+|W^3W^3|)$$

$$+\frac{1}{e^{\beta (B_mJ)}+2+2e^{\beta B_m}}(|W^5W^5|+|W^6W^6|)$$

$$+\frac{1}{1+2e^{\beta (B_mJ)}+2e^{\beta (2B_mJ)}}|111111|$$

$$\{\begin{array}{cc}\frac{1}{2}(|W^5W^5|+|W^6W^6|)\hfill & \mathrm{if}0<B_m<J,\hfill \\ \frac{1}{3}(|W^5W^5|+|W^6W^6|+|111111|)\hfill & \mathrm{if}B_m=J,\hfill \\ |111111|\hfill & \mathrm{if}J<B_m2J\hfill \end{array}$$
(7)
with $`Z=e^{\frac{3}{2}\beta B_m}+2e^{\frac{1}{2}\beta (B_m+2J)}+2e^{\frac{1}{2}\beta (B_m2J)}`$. Eq.(4), in the zero temperature limit, reduces to

$$C(\rho _{AB})=\frac{2}{3(2+e^{\beta (B_mJ)})}\{\begin{array}{cc}\frac{1}{3}\hfill & \mathrm{if}0<B_m<J,\hfill \\ \frac{2}{9}\hfill & \mathrm{if}B_m=J,\hfill \\ 0\hfill & \mathrm{if}J<B_m2J.\hfill \end{array}$$
(8)
## Case III: $`2J<B_m<4J`$

$$\chi _{ABC}=\frac{1}{Z}\left[e^{\frac{1}{2}\beta (B_m+2J)}(|W^5W^5|+|W^6W^6|)+e^{\frac{3}{2}\beta B_m}|111111|\right]|111111|$$
(9)
with $`Z=e^{\frac{3}{2}\beta B_m}+2e^{\frac{1}{2}\beta (B_m+2J)}`$.

## Case IV: $`4JB_m`$

$$\chi _{ABC}=\frac{1}{Z}[e^{\frac{1}{2}\beta (B_m4J)}|W^4W^4|+e^{\frac{1}{2}\beta (B_m+2J)}(|W^5W^5|$$

$$+|W^6W^6|)+e^{\frac{3}{2}\beta B_m}|111111|]|111111|$$
(10)
with $`Z=e^{\frac{3}{2}\beta B_m}+e^{\frac{1}{2}\beta (B_m4J)}+2e^{\frac{1}{2}\beta (B_m+2J)}`$. Clearly, the concurrence in III and IV are both zero. $`B_m=J`$ therefore marks the point of quantum phase transition from an entangled phase to an unentangled one.

It is obvious from Eq.(6) that in Case I, the thermal concurrence remains zero even at nonzero temperatures Wang2 . However, for $`B_m>J`$, unequal mixing of entangled and unentangled states in the spectra of the three-qubit Heisenberg XX ring results in nonzero thermal concurrence at $`T>0`$. The thermal entanglement would, in general, decrease in quantity as $`\chi _{ABC}`$ approaches the maximally mixed state $`\frac{1}{8}I_{ABC}`$ in the limit of infinite temperature. It is thus an interesting problem to determine the critical temperatures $`T_1`$ beyond which the amount of thermal entanglement becomes zero. From Eq.(4), $`T_1`$ clearly depends on the external magnetic field $`B_m`$, in contrast to the two-qubit XX Heisenberg chain Wang1 . These $`T_1`$’s can be obtained by numerically solving

$$2(e^{\beta J}e^{2\beta J})\mathrm{cosh}\frac{1}{2}\beta B_m\sqrt{(e^{2\beta J}+2e^{\beta J})^2+6(e^{2\beta J}+2e^{\beta J})\mathrm{cosh}\beta B_m+9}=0$$
(11)
(see Table I). We note that $`T_1`$ increases with increasing $`B_m`$. When $`B_m`$ is large enough, Eq.(11) reduces to

$$x^66x^52x^33x^2+1=0,$$
(12)
with $`xe^{\beta J}`$, which can be numerically solved to yield $`T_1^{}0.554641J`$. Consequently, as long as $`T<T_1^{}`$, the thermal concurrence would be nonzero, albeit very small.

For $`J<0`$, Eq.(2) reduces, in the zero temperature limit, i.e., $`\beta \mathrm{}`$, to the following four possibilities.

## Case V: $`B_m=0`$

$$\chi _{ABC}=\frac{1}{Z}\left[|000000|+e^{2\beta J}(|W^1W^1|+|W^4W^4|)+|111111|\right]$$

$$\frac{1}{2}(|W^1W^1|+|W^4W^4|)$$
(13)
with $`Z=2+2e^{2\beta J}`$. From Eq.(4), the above equally-weighted mixture has

$$C(\rho _{AB})=\mathrm{max}\{\frac{e^{2\beta J}(3+4e^{\beta J})}{3(1+e^{2\beta J}+2e^{\beta J})},0\}\frac{1}{3}.$$
(14)
## Case VI: $`0<B_m<2J`$

$$\chi _{ABC}=\frac{1}{Z}\left[e^{\frac{1}{2}\beta (B_m+4J)}|W^1W^1|+e^{\frac{1}{2}\beta (B_m4J)}|W^4W^4|+e^{\frac{3}{2}\beta B_m}|111111|\right]|W^4W^4|$$
(15)
with $`Z=e^{\frac{3}{2}\beta B_m}+e^{\frac{1}{2}\beta (B_m4J)}+e^{\frac{1}{2}\beta (B_m+4J)}`$. In the zero temperature limit, Eq.(4) reduces to

$$C(\rho _{AB})=\mathrm{max}\{\frac{2}{3(1+e^{\beta (B_m+2J)})},0\}\frac{2}{3},$$
(16)
a signature of $`|W^4`$.

## Case VII: $`2JB_m4J`$

$$\chi _{ABC}=\frac{1}{Z}[e^{\frac{1}{2}\beta (B_m+4J)}|W^1W^1|+e^{\frac{1}{2}\beta (B_m4J)}|W^4W^4|$$

$$+e^{\frac{1}{2}\beta (B_m+2J)}(|W^5W^5|+|W^6W^6|)+e^{\frac{3}{2}\beta B_m}|111111|]$$

$$\{\begin{array}{cc}\frac{1}{2}(|W^4W^4|+|111111|)\hfill & \mathrm{if}B_m=2J,\hfill \\ |111111|\hfill & \mathrm{otherwise}\hfill \end{array}$$
(17)
with $`Z=e^{\frac{3}{2}\beta B_m}+2e^{\frac{1}{2}\beta (B_m+2J)}+e^{\frac{1}{2}\beta (B_m4J)}+e^{\frac{1}{2}\beta (B_m+4J)}`$. It follows from Eq.(16) and Eq.(17) that

$$C(\rho _{AB})=\{\begin{array}{cc}\frac{1}{3}\hfill & \mathrm{if}B_m=2J,\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}$$

## Case VIII: $`4J<B_m`$

$$\chi _{ABC}=\frac{1}{Z}[e^{\frac{1}{2}\beta (B_m4J)}|W^4W^4|+e^{\frac{1}{2}\beta (B_m+2J)}(|W^5W^5|$$

$$+|W^6W^6|)+e^{\frac{3}{2}\beta B_m}|111111|]|111111|$$
(18)
with $`Z=e^{\frac{3}{2}\beta B_m}+2e^{\frac{1}{2}\beta (B_m+2J)}+e^{\frac{1}{2}\beta (B_m4J)}`$. Consequently, the concurrence is zero in this case. So, $`B_m=2J`$ marks the point of quantum phase transition in the ferromagnetic ring.

The critical temperatures $`T_1`$ for the ferromagnetic ring again depends on $`B_m`$ and can be obtained by numerically solving

$$2(e^{2\beta J}e^{\beta J})\mathrm{cosh}\frac{1}{2}\beta B_m\sqrt{(e^{2\beta J}+2e^{\beta J})^2+6(e^{2\beta J}+2e^{\beta J})\mathrm{cosh}\beta B_m+9}=0$$
(19)
(see Table II). In particular, when $`B_m=0`$, Eq.(19) reduces to

$$e^{2\beta J}(3+4e^{\beta J})=0$$
(20)
which yields $`T_11.27136J`$ Wang2 . We note that $`T_1`$ similarly increases with increasing $`B_m`$. For large enough $`B_m`$, Eq.(19) reduces to

$$y^63y^42y^36y+1=0,$$
(21)
with $`ye^{\beta J}`$, which can be numerically solved, giving $`T_1^{}1.32639J`$. So, for large $`B_m`$, the thermal concurrence would be small but nonzero as long as $`T<T_1^{}`$.

For both $`J>0`$ and $`J<0`$, $`T_1`$ increases with increasing $`B_m`$ up to $`T_1^{}`$ and $`T_1^{}`$ respectively, as long as $`B_m`$ is not infinitely large, in which case $`\chi _{ABC}|111`$ and the thermal concurrence would become zero. Physically, one could understand this phenomenon by looking at the thermal density operator, Eq.(2). Recall that while an equally weighted mixture of $`|W^2`$, $`|W^3`$, $`|W^5`$ and $`|W^6`$ has zero concurrence, Eq.(5) and Eq.(6); an equally weighted mixture of $`|W^5`$ and $`|W^6`$, Eq.(7) and Eq.(8), and $`|W^4`$, Eq.(15) and Eq.(16), have nonzero concurrence. At nonzero temperatures, increasing $`B_m`$ creates a diminishing proportion of $`|W^2`$ and $`|W^3`$, but an increasing proportion of $`|W^4`$, $`|W^5`$, $`|W^6`$ (entangled states), and of course $`|111`$ (unentangled state). The small but nonzero proportion of entangled states contributes to the nonzero thermal concurrence.

Now we describe the quantum teleportation protocol $`P_1`$ using the above three qubit mixed state $`\chi _{ABC}`$ as a resource. It involves a sender, Alice (at site $`A`$), and two receivers, Bob (at $`B`$) and Cindy (at $`C`$). Alice is in possession of two two-level quantum systems, the input system $`S`$, and another system $`A`$ entangled with both a third two-level target system $`B`$ in Bob’s possession, and a fourth two-level target system $`C`$ in Cindy’s possession (i.e. a three-particle entangled state). Here, we label the entangled systems by the site indices. Initially the composite system $`SABC`$ is prepared in a state with density operator

$$\sigma _{SABC}^{total}=\pi _S\chi _{ABC}$$

where

$$\pi _S=|\psi _S\psi |,|\psi _S=\mathrm{cos}\frac{\theta }{2}|0_S+e^{i\varphi }\mathrm{sin}\frac{\theta }{2}|1_S,$$
(22)
$`0\theta \pi ,0\varphi 2\pi `$ are the polar and azimuthal angles respectively, and $`\chi _{ABC}`$, is as given in Eq.(2). To teleport the input state $`\pi _S`$ to Bob’s target system $`B`$ and Cindy’s target system $`C`$, Alice performs a joint Bell basis measurement on systems $`S`$ and $`A`$, described by operators $`\mathrm{\Pi }_{SA}^jI_{BC}`$, $`I_{BC}`$ is the identity operator on the composite subsystem $`BC`$, $`j`$ labels the outcome of the measurement,

$$\mathrm{\Pi }_{SA}^1=|\mathrm{\Phi }^+_{SA}\mathrm{\Phi }^+|,\mathrm{\Pi }_{SA}^2=|\mathrm{\Phi }^{}_{SA}\mathrm{\Phi }^{}|,\mathrm{\Pi }_{SA}^3=|\mathrm{\Psi }^+_{SA}\mathrm{\Psi }^+|,\mathrm{\Pi }_{SA}^4=|\mathrm{\Psi }^{}_{SA}\mathrm{\Psi }^{}|,$$
(23)
where

$$|\mathrm{\Phi }^\pm _{SA}=\frac{1}{\sqrt{2}}(|00_{SA}\pm |11_{SA}),$$

$$|\mathrm{\Psi }^\pm _{SA}=\frac{1}{\sqrt{2}}(|01_{SA}\pm |10_{SA})$$

are the Bell states. If Alice’s measurement has outcome $`j`$, she broadcasts her measurement result (two-bit) to Bob and Cindy via a classical channel. The joint state of Bob’s target system $`B`$ and Cindy’s target system $`C`$ conditioned on Alice’s measurement result $`j`$ is given by

$$\rho _{BC}^j=\frac{1}{p_j}\mathrm{tr}_{SA}[(\mathrm{\Pi }_{SA}^jI_{BC})(\pi _S\chi _{ABC})],$$
(24)
where

$$p_j=\mathrm{tr}_{SABC}[(\mathrm{\Pi }_{SA}^jI_{BC})(\pi _S\chi _{ABC})].$$
(25)
Substituting Eq.(2), Eq.(22), and Eq.(23) into Eq.(25) yields

$$p_1=p_2=\frac{1}{12Z}e^{2\beta (B_m+J)}(f+g\mathrm{cos}\theta ),$$

$$p_3=p_4=\frac{1}{12Z}e^{2\beta (B_m+J)}(fg\mathrm{cos}\theta ),$$
(26)
where

$$f=3e^{\frac{3}{2}\beta B_m}(1+e^{\beta B_m})(1+2e^{3\beta J})+3e^{\frac{1}{2}\beta B_m}(1+e^{3\beta B_m})e^{2\beta J},$$

$$g=e^{\frac{3}{2}\beta B_m}(1e^{\beta B_m})(1+2e^{3\beta J})+3e^{\frac{1}{2}\beta B_m}(1e^{3\beta B_m})e^{2\beta J}.$$

For Bob and Cindy to successfully complete the teleportation protocol, they perform $`j`$-dependent unitary operations $`U_B^j=U_C^j=U^j`$ on systems $`B`$ and $`C`$ respectively $`(\rho _B^j=\mathrm{tr}_C\rho _{BC}^j=\mathrm{tr}_B\rho _{BC}^j=\rho _C^j=\rho ^j)`$ such that

$$\tau _B^j=\tau _C^j=\tau ^j=U^j\rho ^jU^j,$$
(27)
where $`U^j`$ could either be the identity matrix or one of the Pauli matrices (see Table III). The success of the teleportation scheme can be measured by the fidelity Jozsa  between the input state $`\pi _{in}`$ and the output state $`\tau _{out}^j`$, averaged over all possible Alice’s measurement outcomes $`j`$ and over an isotropic distribution of input states $`\pi _{in}`$:

$$F=\frac{1}{4\pi }_0^\pi _0^{2\pi }\mathrm{sin}\theta d\theta d\varphi \underset{j=1}{\overset{4}{}}p_jF^j$$
(28)
where

$$F^j\mathrm{tr}(\tau _{out}^j\pi _{in}).$$
(29)
It follows from Eq.(22) and results from Eq.(24) that

$$F^1=F^2=\frac{h_1+h_2\mathrm{cos}2\theta }{4(f+g\mathrm{cos}\theta )},$$

$$F^3=F^4=\frac{h_1+h_2\mathrm{cos}2\theta }{4(fg\mathrm{cos}\theta )},$$
(30)
where

$$h_1=3e^{\frac{3}{2}\beta B_m}(1+e^{\beta B_m})(3+4e^{3\beta J})+3e^{\frac{1}{2}\beta B_m}(1+e^{3\beta B_m})e^{2\beta J},$$

$$h_2=e^{\frac{3}{2}\beta B_m}(1+e^{\beta B_m})(14e^{3\beta J})3e^{\frac{1}{2}\beta B_m}(1+e^{3\beta B_m})e^{2\beta J}.$$

Substituting Eq.(26) and Eq.(30) into Eq.(28) gives

$$F=\frac{1}{3}+\frac{2}{9}\frac{(2+e^{3\beta J})\mathrm{cosh}\frac{1}{2}\beta B_m}{(1+2e^{3\beta J})\mathrm{cosh}\frac{1}{2}\beta B_m+e^{2\beta J}\mathrm{cosh}\frac{3}{2}\beta B_m}$$
(31)
In order to transmit $`\pi _{in}`$ with fidelity better than any classical communication protocol, we require $`F`$ to be strictly greater than $`\frac{2}{3}`$. In other words, we require

$$\frac{1}{3}(e^{2\beta J}4e^{\beta J})>\frac{\mathrm{cosh}\frac{3}{2}\beta B_m}{\mathrm{cosh}\frac{1}{2}\beta B_m}$$
(32)
and hence $`J<0`$. That is, the nonzero thermal entanglement for $`J>0`$ is “not suitable” as a resource for teleportation via $`P_1`$.

In the zero temperature limit, Eq.(31) reduces to

$$F=\frac{1}{3}+\frac{4}{9(1+e^{\beta (B_m+2J)})}\{\begin{array}{cc}\frac{7}{9}\hfill & \mathrm{if}0B_m<2J,\hfill \\ \frac{5}{9}<\frac{2}{3}\hfill & \mathrm{if}B_m=2J,\hfill \\ \frac{1}{3}<\frac{2}{3}\hfill & \mathrm{if}2J<B_m.\hfill \end{array}$$
(33)
So, in spite of the fact that in Case V, the concurrence is only $`\frac{1}{3}`$, the equally weighted mixture in Eq.(13) is able to yield $`F=\frac{7}{9}`$. Comparing this with the equally weighted mixture of $`|W^5`$ and $`|W^6`$ in Eq.(7), which has concurrence also equal to $`\frac{1}{3}`$, but cannot yield $`F>\frac{2}{3}`$, certainly illustrates a fundamental difference between the entangled mixed states not reflected by the concurrence. For $`0<B_m<2J`$, $`F=\frac{7}{9}`$ is a clear signature of $`|W^4`$ (see Ref.Yeo2 ). At $`B_m=2J`$, the mixing of $`|W^4`$ with an equal proportion of $`|111`$ deteriorates the “quality” of the entanglement so much that $`F`$ is now less than $`\frac{2}{3}`$. We note that $`B_m=2J`$ marks the point of “transition” from $`F>\frac{2}{3}`$ to $`F\frac{2}{3}`$. This coincides with the point of quantum phase transition in the ferromagnetic ring.

For nonzero temperatures it is again an interesting problem to determine the critical temperatures $`T_2`$ beyond which $`F\frac{2}{3}`$. From Eq.(32), $`T_2`$ is clearly dependent on the magnetic field $`B_m`$, as in Ref.Yeo1 . They can be obtained by numerically solving

$$(e^{2\beta J}4e^{\beta J})\mathrm{cosh}\frac{1}{2}\beta B_m3\mathrm{cosh}\frac{3}{2}\beta B_m=0.$$
(34)
(see Table II). Interestingly, when $`B_m=0`$, Eq.(34) reduces to Eq.(20) which thus yields $`T_2=T_11.27136J`$. This means that all nonzero thermal entanglement, in this case, is “suitable” as a resource for teleportation via $`P_1`$. The mixing of states here clearly does not have a devastating effect on the quality of the thermal entanglement. Supppose $`T`$ is small enough, and let $`B_m=\eta J`$, $`0<\eta `$, then Eq.(34) yields

$$T_2=\frac{1}{\mathrm{ln}3}(2\eta )J0.$$
(35)
So, $`\eta `$ can at most equal 2, and as $`\eta 2`$, $`T_20`$ consistent with our assumption. As shown in Table II, $`T_2`$ decreases with increasing $`B_m`$ and each $`T_2`$ is strictly less than the corresponding $`T_1`$, which increases asymptotically to $`T_1^{}`$. This means that with increasing $`B_m`$ we have an increasing range of nonzero thermal entanglement which is however not able to yield $`F>\frac{2}{3}`$. Physically, one could attribute the cause of the poor quality of thermal entanglement to the fact that there is now a comparable or greater proportion of unentangled $`|111`$ than the “teleportation grade” $`|W^1`$ and $`|W^4`$.

In conclusion, our “average teleportation fidelity criterion”, Eq.(32) reveals several interesting aspects of quantum entanglement not reflected by concurrence. On the one hand, whereas comparable mixing of entangled states certainly decreases the resulting pairwise concurrence (see Eq.(5) and Eq.(6), Eq.(7) and Eq.(8), Eq.(13) and Eq.(14)), it could result either in low quality states which yield $`F\frac{2}{3}`$, or in high quality states giving $`F>\frac{2}{3}`$. On the other hand, comparable mixing of entangled with unentangled states not only certainly decreases the resulting pairwise concurrence (see Eq.(7) and Eq.(8), Eq.(17)) but definitely degrades the teleportation quality of the entangled mixed state. Furthermore, the teleportation quality of the entangled mixed state is more sensitive to the degree of mixing than its concurrence. As a result, we could have a more entangled thermal state not giving a better average fidelity than a less entangled one. Since entanglement is such an important resource in quantum information, it is very important to have a more fundamental understanding of these aspects of quantum entanglememt.

The author thanks Yuri Suhov and Andrew Skeen for useful discussions. This publication is an output from project activity funded by The Cambridge MIT Institute Limited (“CMI”). CMI is funded in part by the United Kingdom Government. The activity was carried out for CMI by the University of Cambridge and Massachusetts Institute of Technology. CMI can accept no responsibility for any information provided or views expressed.


format: text

text:
# On the law of motion in Special Relativity

## Abstract

Newton’s law of motion for a particle of mass $`m`$ subject to a force $`𝒇`$ acting at time $`t`$ may be formulated either as

$$𝒇=d/dt(m𝒖(t))$$

or, since $`m=m_0`$ is a constant, as

$$𝒇=m𝒂(t)$$

where $`𝒖(t)`$ and $`𝒂(t)`$ are the velocity and the acceleration, respectively, of the particle at time $`t`$ relative to an inertial frame $`S,`$ ‘the laboratory’. This law may be interpreted in either of two ways:

1. The force $`𝒇`$ acting on the particle at time $`t`$ during an infinitesimal time $`\delta t`$ imparts to the laboratory a boost $`\delta 𝒖=(1/m)𝒇\delta t,`$ while the particle maintains the velocity $`𝒖(t)`$ relative to the new frame $`S^{}`$.
2. The force $`𝒇`$ acting on the particle at time $`t`$ during an infinitesimal time $`\delta t`$ imparts to the particle a boost $`\delta 𝒘=(1/m)𝒇\delta t`$ relative to its proper frame $`S_0`$ which moves with velocity $`𝒖(t)`$ relative to $`S`$.

We show that the relativistic law of motion admits both interpretations, the first of which is in fact equivalent to the law of motion. As a consequence, we show that the relativistic law of motion may also be reformulated as

$$𝒇=m𝒂$$

in analogy with Newton’s law, but with a relativistic mass and a relativistic acceleration defined in terms of the relativistic addition law of velocities, rather than ordinary mass and ordinary vectorial addition of velocities that lead to the classical acceleration and to Newton’s law.

## Introduction

It is well-known that the Special Theory of Relativity is based on two postulates: the principle of relativity and the constancy of the speed of light. From these postulates, one can deduce the Lorentz transformations that connect the space-time coordinates of a particular event in two inertial frames. From these transformations, in turn, lenght contraction, time dilation, the ‘addition’ formula for velocities, etc., follow straightforwardly (see e.g , , , ). Thus, these postulates are sufficient for the development of the theory as far as kinematics is concerned.

On the other hand, this is not so for dynamics. In the earlier papers, the law of motion was obtained from electrodynamical considerations. After a tortuous path, Planck and independently Tolman, finally arrived at the familiar formulation of today connecting force, mass and velocity at a given time:

$$\{\begin{array}{c}𝒇=d/dt(m𝒖(t))\\ m=m_0\gamma (u)\\ \gamma (u)=1/\sqrt{(1u^2/c^2)}\end{array}$$

The time $`t`$, the mass $`m`$, the force $`𝒇`$ and the velocity $`u(t)`$ are relative to an inertial frame $`S`$. (See , ,  and  for historical references).

In his lecture on ‘Space and Time’ (see ), delivered in 1908 and published a year later, Minkowski introduced an invariant reformulation of the usual law of motion:

$$𝑭=m_0𝑨$$

where $`𝑭`$ is the 4-force acting on the particle and $`𝑨`$ the 4-acceleration.

It should be noticed, however, that $`𝑭`$ is a 4-force and thus differs from the 3-force $`𝒇`$. Even the spatial component $`(F_1,F_2,F_3)`$ of Minkowski 4-force differs from the previously introduced 3-forces. In fact

$$(F_1,F_2,F_3)=\gamma (u)𝒇$$

The corresponding spatial component of the 4-acceleration is therefore

$$(A_1,A_2,A_3)=\gamma (u)d/dt(\gamma (u)𝒖)$$

The time components of the Minkowski force and the acceleration are given by $`F_4=(\gamma (u)/c)𝒇.𝒖`$ and $`A_4=c\gamma (u)d/dt(\gamma (u)).`$ The equation $`F_4=m_0A_4`$ follows from the corresponding equation for the spatial component and is in fact equivalent to $`𝒇.d𝒓=c^2dm,`$ as shown by an easy calculation.

The arguments generally given in favor of the usual law of motion are (see ):

1. it reduces to the law of Newton $`𝒇=m𝒂`$ in the classical limit
2. it leads to the conservation of momentum in simple collisions, provided that we assume the law of equality of action and reaction at contact
3. it leads to $`\delta E=c^2\delta m,`$ the infinitesimal version of Einstein’s famous law $`E=mc^2`$ by adopting the classical definition of work as $`force\times distance,`$ as we will see later on.
4. it is consistent with the well-established Lorentz law of force in electrodynamics

These are strong arguments not only in favor of the usual law of motion, but (iii) and (iv) provide good reasons to consider also the 3-force (along with the Minkowski 4-force ).

The aim of this note is to show that this law of motion admits the two interpretations of Newton’s law stated in the Abstract. To formulate them in a concise manner, we shall use the diagram

as a shorthand to describe the following situation: the particle $`P`$ moves with velocity $`𝒖^{}`$ relative to $`S^{},`$ which itself moves with velocity $`𝒗`$ relative to $`S`$ and $`𝒖`$ is the velocity of the particle relative to $`S.`$ It will be particularly handy when we deal with relativity and Lorentz transformations.

To simplify the notation, we shall use either $`_t`$ or $`\dot{(\mathrm{})}`$ for $`d/dt`$. i.e., the derivative with respect to time.

## Newtonian Dynamics

Newton’s law of motion for a particle of mass $`m`$ subject to a force $`𝒇`$ acting at time $`t`$ is

$$𝒇=_t(m𝒖(t))$$

where $`𝒖(t)`$ is the velocity of the particle at time $`t`$ relative to an inertial system $`S,`$ to be referred to as ‘the laboratory’.

This law can be rewritten as

$$𝒖(t+\delta t)=𝒖(t)+(1/m)𝒇\delta t$$

by letting $`\delta t`$ an infinitesimal time increment and using

$$𝒖(t+\delta t)=𝒖(t)+\dot{𝒖}(t)\delta t$$

Now, the point is that this law may be interpreted in either of two equivalent ways:

1. The force $`𝒇`$ acting on the particle at time $`t`$ during an infinitesimal lapse of time $`\delta t`$ imparts to the laboratory $`S`$ , relative to which the particle has velocity $`𝒖`$(t), a boost $`\delta 𝒖=(1/m)𝒇\delta t`$ while the particle maintains the velocity $`𝒖(t)`$ relative to the new frame $`S^{}`$. Thus, the new velocity of the particle relative to $`S`$ at time $`t+\delta t`$ is

$$𝒖(t+\delta t)=(1/m)𝒇\delta t+𝒖(t)$$

Using our diagram,
2. The force $`𝒇`$ acting on the particle at time $`t`$ during an infinitesimal lapse of time $`\delta t`$ imparts to the particle a boost $`\delta 𝒘=(1/m)𝒇\delta t`$ relative to its proper frame $`S_0`$ which moves with velocity $`𝒖(t)`$ relative to $`S`$. Then the new velocity of the particle relative to $`S`$ at time $`t+\delta t`$ is

$$𝒖(t+\delta t)=𝒖(t)+(1/m)𝒇\delta t$$

Diagrammatically,

## Relativistic Dynamics

As we mentioned already, in Special Relativity Newton’s law of motion is replaced by

$$\{\begin{array}{c}𝒇=_t(m𝒖(t))\\ m=m_0\gamma (u)\end{array}$$

From now on, we choose a system of unities such that $`c=1`$. Thus $`\gamma (u)=1/\sqrt{1u^2}`$.

We wish here to point out a suggestive interpretation of the law of motion which does not seem to have been noticed before. In fact, we wish to point out that this law admits precisely the two previous interpretations, provided that we take care to refer all quantities to the corresponding frames, since force, mass and time are frame-dependent, contrary to the classical case, and take into account the fact that relativistic composition of velocities $``$ is not commutative. Unlike the Newtonian case, however, only the first interpretation is equivalent to the law of motion. The second, although a consequence of this law, does not seem to be equivalent to it.

Thus,

1. The force $`𝒇`$ acting on the particle at time $`t`$ during an infinitesimal lapse of time $`\delta t`$ imparts to the laboratory $`S`$ , relative to which the particle has velocity $`𝒖`$(t), a boost $`\delta 𝒖=(1/m)𝒇\delta t`$ (where $`𝒇`$, $`\delta t`$ and m are measured in $`S`$), while the particle maintains the velocity $`𝒖(t)`$ relative to the new frame $`S^{}`$. Thus, the new velocity of the particle relative to $`S`$ at time $`t+\delta t`$ is

$$𝒖(t+\delta t)=(1/m)𝒇\delta t𝒖(t)$$

Using our diagram once again, we obtain
2. The force $`𝒇`$ acting on the particle at time $`t`$ during an infinitesimal lapse of time $`\delta t`$ imparts to the particle a boost $`\delta w_0=(1/m_0)𝒇_0(\delta t)_0`$ relative to $`S_0,`$ its rest frame which moves with velocity $`𝒖(t)`$ relative to $`S`$. (Here $`𝒇`$, $`(\delta t)_0`$ and $`m_0`$ are measured in $`S_0`$). Then the new velocity of the particle relative to $`S`$ at time $`t+\delta t`$ (measured in $`S`$) is

$$𝒖(t+\delta t)=𝒖(t)(1/m_0)𝒇_0(\delta t)_0$$

Diagrammatically,

To prove our claim, we use some basic facts about relativistic kinematics that can be found in the Appendix.

When the velocity $`𝒗`$ is infinitesimal, $`𝒗=\delta 𝒘,`$ we get the infinitesimal Lorentz transformation

$$\begin{array}{ccc}\delta 𝒗𝒖\hfill & =\hfill & (𝒖+\delta 𝒗)/(1+\delta 𝒗.𝒖)\hfill \\ & =\hfill & (𝒖+\delta 𝒗)(1\delta 𝒗.𝒖)\hfill \\ & =\hfill & 𝒖+\delta 𝒗𝒖(𝒖.\delta 𝒗))\hfill \end{array}$$

since $`\gamma (\delta 𝒗)=1`$

Then

$$𝒖._t\gamma u=𝒖.(\dot{\gamma }𝒖+\gamma \dot{𝒖})=\dot{\gamma }u^2+\gamma 𝒖.\dot{𝒖}=\dot{\gamma }u^2+\dot{\gamma }/\gamma ^2=\dot{\gamma }()$$

since $`\dot{\gamma }=\gamma ^3𝒖.\dot{𝒖}.`$

Assuming the law of motion $`𝒇=m_0_t(\gamma 𝒖)`$,

$$\delta w=\delta t𝒇/m=\delta t\gamma ^1_t\gamma 𝒖=\dot{𝒖}\delta t+\delta t(\dot{\gamma }/\gamma )𝒖$$

and

$$𝒖.\delta w=𝒖.\delta t\gamma ^1_t\gamma 𝒖=\delta t\gamma ^1𝒖._t\gamma 𝒖=\delta t(\dot{\gamma }/\gamma )$$

Thus,

$$\begin{array}{ccc}\delta 𝒘𝒖\hfill & =\hfill & 𝒖+\delta 𝒘𝒖(𝒖.\delta w)\hfill \\ & =\hfill & 𝒖+\delta 𝒘𝒖\delta t(\dot{\gamma }/\gamma )\hfill \\ & =\hfill & 𝒖+\dot{𝒖}\delta t\hfill \\ & =\hfill & 𝒖(t+\delta t)\hfill \end{array}$$

We have obtained thus the first interpretation.

Notice that from (\*) we also have

$$𝒖._t(m𝒖)=m_0\dot{\gamma }=\dot{m}$$

so that

$$0=\delta t𝒖.(𝒇_tm𝒖)=\delta E\delta m$$

the infinitesimal version of Einstein’s law. (Recall that $`c=1`$).

To obtain the second interpretation from the relativistic law of motion, we define left and right accelerations by the formulas

$$\{\begin{array}{c}𝒂_l\delta t𝒖(t)=𝒖(t+\delta t)\\ 𝒖(t)𝒂_r\delta t=𝒖(t+\delta t)\end{array}$$

From these equations we obtain unique solutions

$$\{\begin{array}{c}𝒂_l=𝒂+𝒖\gamma ^2(𝒂.𝒖)\hfill \\ 𝒂_r=\gamma 𝒂+𝒖(1/𝒖^2)(\gamma (\gamma 1)𝒂.𝒖)\hfill \end{array}$$

Furthermore, it is easily checked that

$$\{\begin{array}{c}𝒂_l=(1/\gamma )[𝒂_r+𝒖(1/𝒖^2)𝒂_r.𝒖(\gamma 1)]\hfill \\ 𝒂_r=\gamma (𝒂_l+𝒖[(1/𝒖)^2𝒂_l.𝒖(\gamma 1)𝒂_l.𝒖\gamma ])\hfill \end{array}$$

Now, $`(\delta t)_0=(1/\gamma )\delta t`$ (time dilation) and thus,

$$(1/\gamma m_0)𝒇_0=𝒂_riff𝒇_0=\gamma \{𝒇+𝒖[(1/𝒖)^2(𝒇.𝒖(\gamma 1))𝒇.𝒖\gamma ]\}$$

But from the law of motion $`𝒇=_t(m𝒖)`$ (and $`m=m_0\gamma `$), it follows easily that $`𝒇=m𝒂_l`$ which in turn, implies the right hand side. But the right hand side is true, since it is precisely the Lorentz transformation of the force, as deduced from the law of motion (see the Appendix)

To show that, conversely, the first interpretation (together with $`m=m_0\gamma `$) implies the relativistic law of motion, notice that from the first interpretation, $`(1/m)𝒇\delta t=𝒂_l\delta t,`$ which in turn implies $`𝒇=m𝒂_l,`$ a formula which is obviously equivalent (by taking derivatives) to $`𝒇=_t(m𝒖(t),`$ provided that $`m=m_0\gamma .`$

Thus, we may reformulate the relativistic law of motion as

$$𝒇=m𝒂_l$$

where $`m=m_0\gamma `$ is the relativistic mass of the particle and $`𝒂_l`$ is the relativistic ”left acceleration”.

The second interpretation, $`(1/\gamma m_0)𝒇_0=a_r,`$ can be written in a way that is analogous to Newton’s law, namely

$$𝒇_0=m_0\delta 𝒘_0/(\delta t)_0$$

(since $`\delta 𝒘_0=1/m_0𝒇_0(\delta t)_0`$), but with the relativistic ”acceleration” $`\delta 𝒘_0/(\delta t)_0=\gamma 𝒂_r`$. This expression has the following physical interpretation: assume that an observer attached to the frame $`S_0`$ moving with the particle let a test body ‘fall’ freely relative to $`S`$. Then $`\delta 𝒘_0/(\delta t)_0`$ is the acceleration at the moment of the take-off, as measured by this observer in his frame (i.e. $`S_0`$). In classical mechanics, this acceleration is implicit in the studies of Huygens on centrifugal force (De vi centrifuga). In fact, Huygens imagines a man attached to a turning wheel and holding a thread tied to a ball of lead in his hand. The thread is suddenly cut and Huygens studies the motion of the ball at the instant when the thread is cut. This ”take-off” acceleration plays an important role in some historico-critical studies such as those of the ‘Ecole du fil’ of of F.Reech and J.Andrade (see ).

## Appendix:Lorentz transformations and composition of velocities

We recall the Lorentz transformation that connects the position of an event in two inertial frames $`S`$ and $`S^{}`$ such that $`S^{}`$ (or rather its origin) moves with uniform velocity $`𝒗`$ with respect to $`S.`$ Indeed, if $`(𝒓,t)`$ describes the event in $`S,`$ and $`(𝒓^{},t^{})`$ describes the same event in $`S^{},`$ then (by mapping these quantities in an independent Euclidean 3-space):

$$\{\begin{array}{c}𝒓^{}=𝒓+𝒗\gamma (v)\{(1/c^2)\gamma (v)/(1+\gamma (v))𝒗.𝒓t\}\\ t^{}=\gamma (v)(t𝒗.𝒓/c^2)\end{array}$$

Similarly, if a particle moves with velocity $`𝒖`$ with respect to $`S`$ and velocity $`𝒖^{}`$ with respect to $`S^{},`$ then

$$𝒖^{}=\{𝒖/\gamma (v)+𝒗[(1/c^2)\gamma (v)/(1+\gamma (v))𝒖.𝒗1]\}/(1𝒖.𝒗/c^2)$$

Finally, if the particle is subject to a force $`𝒇`$ relative to $`S,`$ then relative to $`S^{}`$ the force is given by

$$𝒇^{}=\{𝒇/\gamma (v)+𝒗[(1/c^2)\gamma (v)/(1+\gamma (v))𝒇.𝒗𝒇.𝒖/c^2]\}/(1𝒖.𝒗/c^2)$$

These formulas may be found in  (pages 23, 40 and 97), although the expression $`(11/\gamma (v))/v^2`$ has been replaced by $`(1/c^2)\gamma (v)/(1+\gamma (v))`$ used by Ungar . This last formula makes sense for every $`𝒗`$ such that $`v<c.`$

The Lorentz transformation of velocities may be expressed in terms of composition or addition of velocities.

In fact, let $`V=\{\stackrel{}{v}\text{R}^3|v^2<c^2\}.`$ Define an operation on the open domain $`V\times V`$ of $`\text{R}^3\times \text{R}^3`$ with values in $`V`$ by the formula

$$𝒗𝒖=\{𝒖/\gamma (v)+𝒗[(1/v^2)(11/\gamma (v))𝒖.𝒗+1]\}/(1+𝒖.𝒗/c^2)$$

Although this operation is neither commutative nor associative, it is gyrocommutative and gyroassociative. (See  for these notions as well as for further properties of this operation).

In terms of this operation, we can express the Lorentz transformation of velocities of the beginning of this Appendix as

$$𝒖^{}=(𝒗)𝒖$$

or, equivalently, as

$$𝒖=𝒗𝒖^{}$$

Notice that even if $`𝒖`$ and $`𝒗`$ are velocities, $`𝒗𝒖`$ is not interpretable as a velocity. To make contact with actual composition of velocities we need to devise physical set-ups to realize each vector as a velocity (relative to a suitable frame) and the operation as actual composition of velocities according to our diagram

## Acknowledgments

The first author owes a great debt of gratitude to Nicanor Parra. He formulated clearly the problem of the status of the 3-force and the law of motion in the Theory of Special Relativity and told him about the take-off acceleration in classical mechanics, suggesting that it could be used in Special Relativity. He had countless conversations on and off on ‘natural and supernatural’ mechanics, since he followed Parra’s course in 1957 at the Instituto Pedagógico of the Universidad de Chile in Santiago. His encouragement is greatly appreciated. He is also in debt to Jorge Krausse. Besides discussions of a general nature, he helped him to find his way in the literature on Relativity Theory.


text:
# Event-Driven Molecular Dynamics in Parallel

## I Introduction

Event-driven molecular dynamics is an effective algorithm for the simulation of many-component systems, which evolve independently, except for discrete asynchronous instantaneous interactions. As an example we will discuss a system consisting of $`N`$ hard spheres in a box with periodic boundary conditions, but the algorithm can be extended to particles with different shapes and interacting via any piecewise constant potential or to completely different problems as well lubachevsky00 .

Event-driven molecular dynamics processes a series of events asynchronously one after another. A straight-forward but simplistic approach alder59  updates all particles at each event; many advanced versions are available in text-books allen87 ; rapaport95 . We base our algorithm on the sophisticated algorithm as presented by lubachevsky91 , which updates only those particles involved in an event. It has been succesfully applied to many different problems, among them granular gases luding99 ; luding02c , polymer chains smith97 , tethered membranes luding02c ; sm:memb , protein folding dokholyan98 , and battlefield models nicol88 .

In many physical systems the duration of the interaction of the components, e. g. the collision of two particles, is negligible compared to the time between these interactions. The simulation of such systems with traditional time-driven molecular dynamics is highly inefficient. Instead, it is straight-forward to consider the interactions as instanteneous events and solve the problem with an event-driven algorithm.

One reason why event-driven molecular dynamics has not become as popular as conventional time-driven molecular dynamics is the fact that parallelization is extremely complicated. The paradoxical task is to algorithmically parallelize physically non-parallel dynamics. Nevertheless some ideas and basic considerations about the parallelization have been proposed in lubachevsky92 . They are especially suited for shared memory computers, but can be transferred to distributed memory architectures as well marin97 . Apart from those ideas no full and general implementation of a parallelized algorithm including load-balancing has been published yet, to our knowledge. In this paper we present an algorithm, which is based on those ideas, but is enhanced and completed at several points. It can be implemented with generic tools such as MPI, and therefore it is suited for shared and distributed memory architectures alike.

In section II we explain the details of the implementation of event-driven molecular dynamics. The parallelization of the algorithm is presented in section III and a summary is given in section IV.

## II Event-Driven Molecular Dynamics

Section II.1 gives an outline of the main routine of the algorithm, which is composed of 4 steps (see Fig. 1). Then the most important data structures are introduced in section II.2 and hereby step 1 and 4 are treated. Step 2 and 3 are discussed in sections II.3 and II.4, respectively. Finally, section II.5 deals with performance issues.

### II.1 Outline

Event-driven molecular dynamics processes a series of discrete events. In a system of hard spheres events typically refer to instantaneous collisions, involving exactly two particles. Only these two particles are processed; the state of the other particles remains unchanged. So the state information of most particles refers not to the current time, but to a different point of time in the past.

Event processing includes a state update for the concerned particles (see section II.3) and the calculation of the next future events for those particles (see section II.4). An outline of the algorithm can be seen in Fig. 1.

The serial algorithm only adds to this event processing loop an initialization step at the beginning and periodic data output. The parallel algorithm needs several additional routines, which are described in section III.

### II.2 Data Organization

The algorithm maintains two data fields, particle states and event list, which contain exactly one entry per particle (see Fig. 2). The former refers to the past of a particle, the latter to its future ($`i:\text{state time }t_0(i)\text{current time}\text{event time }t_{\mathrm{ev}}(i)`$).

A particle state consists of the physical state of a particle, such as position, velocity, etc. immediately after the most recently processed event of that particle, the point of time of that event, an event counter (see section II.3), a cell number (see section II.4), and a global particle number, which is needed to identify particles on different processes (see section III).

The event list associates with every particle an event in the future. The data units of the event list consist of event time, event type, event partner, and a copy of the event counter of the partner.

All the events are scheduled in a priority queue which is usually implemented as an implicit heap tree lubachevsky91 ; lubachevsky92 ; knuth68 , but other data structures with similar characteristics are possible, too marin95 ; knuth68 . A heap tree is a binary tree with the following ordering rule: Each parental node has higher priority than its children; the children themselves are not ordered. So the root node always has highest priority, i. e. the earliest event time. To get the next event from the priority queue takes computational costs of $`𝒪\left(1\right)`$. Insertion of a new event in the priority queue is done with costs of $`𝒪\left(\mathrm{log}N\right)`$.

This data organization is more efficient compared to the algorithm in lubachevsky91 ; lubachevsky92 , since no double buffering is needed here, however, the basic ideas of lubachevsky91 ; lubachevsky92  are still valid.

### II.3 State Update

When a collision between two particles is processed, the states of these particles are updated from a point of time in the past to the current simulation time. First, the positions and velocities of the particles immediately before the collision can be derived by inserting the old state in the equation of motion. As a result of this first step the particles are touching each other. Then the interaction between the particles takes place and yields new particle velocities. Now, the event counter is increased, the state time is updated to the current simulation time, and the values of the positions and velocities immediately after the collision are stored in the state array.

A typical collision rule for hard spheres luding01d ; luding02b  looks like

$`𝐯_{}^{}{}_{1/2}{}^{}`$ $`=𝐯_{1/2}{\displaystyle \frac{1+r}{2}}\left(\widehat{𝐤}(𝐯_1𝐯_2)\right)\widehat{𝐤},`$
where primes indicate the velocities $`𝐯`$ after the collision, $`\widehat{𝐤}`$ is a unit vector pointing along the line of centers from particle 1 to particle 2, and $`r`$ denotes the restitution coefficient. For the performance tests in the subsequent chapters we have used the simplest case without dissipation ($`r=1`$). 111Soft spheres can be approximated via particles consisting of concentric shells with piecewise constant interaction potentials.

If the event partner has undergone collisions with other particles between the scheduling and the processing of the event, it becomes invalid. Validity of event partners can be checked by comparing the event counter of the partner to the copy in the event list. If they do not match, the partner has collided with other particles in the meantime, and the scheduled collision is no longer valid. Then event processing only refers to one particle and the state update described above only consists of the particle motion; no collision is performed. After that the algorithm continues in the normal way, i. e. the next event for the particle will be calculated and scheduled in the priority queue.

Note that the algorithm in lubachevsky91 ; lubachevsky92  uses a different strategy to detect invalid event partners: The algorithm checks at each collision if an event partner becomes invalid and if so marks this partner. This strategy is less efficient, because sometimes the same partner is invalidated several times. Besides, in the parallel algorithm additional communication between different processes might be necessary.

### II.4 Event Calculation and Linked Cell Structure

When an event has been processed, new events for the particles involved in that event have to be calculated. In simulations of hard spheres this means the calculation of possible future collisions. If the particles move on ballistic trajectories

$`𝐫_i(t)=𝐫_i(t_0)+𝐯_i(t_0)(tt_0)+{\displaystyle \frac{1}{2}}𝐠(tt_0)^2,`$
two particles 1 and 2 will collide at time $`t_{12}`$:

$`t_{12}`$ $`=t_0+\left(𝐫_{12}𝐯_{12}\sqrt{(𝐫_{12}𝐯_{12})^2\left(r_{12}^2(R_1+R_2)^2\right)v_{12}^2}\right)/v_{12}^2,`$
where $`𝐯_{12}=𝐯_2(t_0)𝐯_1(t_0)`$ and $`𝐫_{12}=𝐫_2(t_0)𝐫_1(t_0)`$ are the relative velocities and positions of the particles at time $`t_0`$, and $`R_i`$ are the radii of the particles. If $`t_{12}`$ is imaginary or smaller than $`t_0`$, the particles will not collide.

If the algorithm would have to check for collisions with all other particles, performance would be very poor for large numbers $`N`$ of particles. So we divide simulation space in $`C`$ cells with equal sides. Each particle belongs to the cell in which its center lies. If the cell size is larger than the maximal interaction length of the particles, i. e. the particle diameter in the case of hard spheres, event calculation has to check for possible collisions of two particles only if they belong to the same cell or if their cells are neighbors. (One square cell has 8 neighbors in 2D and a cube has 26 neighbors in 3D.)

However, the algorithm has to keep track of cell changes. So additional events, namely cell changes, come into play. They are treated just in the same way as the collision events; only the collision partner is the boundary of a cell. The difference with a collision event is that only one particle is involved in a cell change, and the velocity of the particle does not change at event time. Cell changes at the boundary of the simulation space require that the position vector jumps to the opposite side of the simulation volume due to the periodic boundary conditions.

### II.5 Optimal Cell Numbers

On average there are $`3^DN/C1`$ particles in the neighborhood of each particle. So in the limit $`CN`$ this number becomes very large and many possible events have to be calculated after each collision. On the other hand, if $`CN`$, then each particle has to cross many cell boundaries between a collision of two particles, and thus more events have to be treated to complete the simulation. These two contributions compete with each other, the first one is proportional to $`(C/N)^1`$ and the second one is proportional to $`(C/N)^{1/D}`$. Fig. 3 (left) shows a broad minimum of the simulation time for low densities at the optimal cell number $`C/N1.5`$ in 2D and $`C/N8`$ in 3D. So the program can choose an optimal cell number before the calculation starts. Note that in the high density limit, especially in 3D, the optimal cell number cannot be reached, because the size of the cells has to be larger than the particle diameter (see Fig. 3 (right)). So, in this case $`C`$ should simply be chosen as big as possible.

## III The Parallel Algorithm

In section III.1 we demonstrate how parallelization can be achieved via domain decomposition and dynamic load-balancing. Then, in section III.2 we point out the difficulties that arise in parallelizing this algorithm and the necessity of state saving and error recovery (see section III.3). In section III.4 the parallel algorithm is explained in detail. Finally, some performance issues are discussed in section III.5.

### III.1 Domain Decomposition and Dynamic Load-Balancing

Parallelization is achieved via domain decomposition. Each cell is affiliated with a process, but this affiliation can change during the simulation if the load of the processes has to be rebalanced.

In order to realize maximal performance, the computational load should be distributed equally among the processes. If inhomogenities exist from the very beginning or emerge during the simulation, dynamic load-balancing becomes important. A typical example is cluster formation in granular gases. luding99 ; luding02

Domain decomposition should thus take into account the following points: Firstly, the load of the processes should be distributed homogeneously. Secondly, in order to minimize process communication, the border area of each process domain should be as small as possible, which implies that the ideal process domains are squares or cubes. Thirdly, a simple and fast function that assigns a cell to a process and vice versa is required.

We meet these demands in the following way: Cell numbers are assigned to cell coordinates in a tree-like structure, see a small 2D-example in Fig. 4. (Note that a realistic example has thousands or millions of cells.) These cell numbers are obtained by interleaving all the bits $`1`$ to $`k`$ of the cell coordinates: $`z[k],y[k],x[k],\mathrm{},z[1],y[1],x[1]`$. The result is a cell number in the range $`0`$ to $`2^{kD}1`$. (For convenience these numbers are increased by one in Fig. 4.) 222An example: The bottom left cell in Fig. 4 with the coordinates 0/3 (bit pattern 000/011) has the number corresponding to the bit pattern 001010, i. e. 10. Incrementation by one yields the cell number 11.

Then a block of $`2^n`$ consecutive cell numbers is affiliated with each process, where $`n`$ can be different for each process. Suppose the cells 1-16 in Fig. 4 should be distributed over 4 processes. If a lot of particles are aggregated in the lower right corner, load-balancing could result in the following layout: Process I is affiliated with cells 1-8, process II with cells 9-12, process III with cells 13-14, and process IV with cells 15-16.

Every now and then the load of the processes has to be checked. A reasonable and simple measure of the load is the number of particles or the number of collisions in the process domain. If a restructuring of the domain decomposition would result in a significantly better load-balance, three processes, respectively, exchange information about their particle and cell data, so that two light-weight neighboring processes can merge their domains and a heavy-weight process can split its domain in two halves. In the example above, if the initially inhomogeneous system becomes homogeneous, this procedure could result in the merging of the domains of processes III and IV and a splitting of the domain of process I. Process II does not participate in this rebalancing, but it should be informed that its neighbor processes have changed, of course. In the end the layout would look like that: Process I is affiliated with cells 1-4, process II with cells 9-12, process III with cells 5-8, and process IV with cells 13-16.

### III.2 Causal Order

A parallel approach to simulate asynchronous events has to make use of the concept of local times. Each process handles the events in its domain and thereby increases the local simulation time. When particles cross domain boundaries, the affected processes communicate with each other and events are inserted into or removed from the event lists of the processes. If the event time of a newly inserted particle is less than the local simulation time, causality is violated, because a collision of that particle with another one, which has already been processed, could have been missed.

In general, there are two strategies to circumvent this problem: In a conservative approach only those processes that are guaranteed not to violate causality are allowed to handle their events, the rest of the processes has to idle. In an optimistic approach the processes do not have to idle. If causality is violated, a complex rollback protocol which corrects the erroneous computations has to be executed. Whether a conservative approach is efficient or not depends highly on the maximal event propagation speed.

Unfortunately, in event-driven molecular dynamics there is no upper limit for the speed of events propagating along a chain of particles lubachevsky92 , even if the particles themselves are moving very slowly. In other words, in the conservative case one process is operating and all others are idling and we are back to the serial algorithm. So we are left with the optimistic strategy and have to undertake the task of implementing a rollback protocol.

### III.3 State Saving and Error Recovery

When a causality error is detected, the simulation is restarted from the latest saved state. So the algorithm has to make a backup copy of the simulation data periodically and has to ensure that there is no latent causality error in this backup. The latter is guaranteed if all processes are synchronized at saving time. This means that only those processes are allowed to continue their computations whose local simulation times have not yet reached synchronization time. Note that there are other operations, like e. g. data output or load-balancing, which require periodic synchronization anyway. Besides, without synchronization the local simulation times tend to drift apart, which makes causality errors very likely korniss01 .

To find the optimal backup interval, we make use of an adaptive strategy. If no causality error turns up between two successive save operations, the interval increases, otherwise it decreases. If other operations trigger synchronization, this occasion is used for state saving, too, of course.

If a causal error occurs, all processes perform a rollback to the saved state. Furthermore a synchronization barrier is scheduled for the time when the error occured. This prevents the same causality error from happening again. Of course, another error could occur at an earlier point of time. Then the simulation would perform another rollback and an earlier resynchronization would be scheduled.

For comparison the algorithm described in lubachevsky92  needs two backup copies of the simulation data. So our strategy reduces memory consumption further by a factor $`3/2`$.

### III.4 Border Zone and Process Communication

The border zone (in lubachevsky92  it is called insulation layer) consists of the cells $`C_{\mathrm{border}}`$ whose neighbors belong to a different process. With the domain decomposition described in section III.1, there is always a “monolayer” of border zone cells at the boundary of a process area.

Each process thus maintains a list of virtual border zone cells which actually belong to its neighboring processes and a list of virtual particles residing in the virtual border zone cells. Those virtual particles can act as partners for real particles. But no events are calculated for them directly and they are not represented in the event list, since they are real particles on another process. The events are already calculated there and will be communicated to the adjacent processes.

However, it is highly inefficient for a process to communicate after every event with its neighbors. So the parallel algorithm is designed in a stepwise manner (see Fig. 5): Each computing step lasts until the next event in the border zone of a process, which is obtained in the following way: An event is associated with each particle. Each event belongs to a cell (or two adjacent cells, if two colliding particles reside in different cells or if a particle changes from one cell to another). The smallest event time in a cell is called the cell time $`t_{\mathrm{cell}}`$. These cell times are stored in an additional priority queue similar to the event list in section II.2. This list contains the scheduled cell times for all local border zone cells and returns the minimum time $`t_{\mathrm{step}}=\mathrm{min}_c(t_{\mathrm{cell}}(c))=t_{\mathrm{cell}}(c^{})`$ which belongs to cell $`c^{}`$.

Now, this time is used as the preliminary length of the calculation step. But firstly, the neighboring cells of $`c^{}`$, if located on other processes are checked. If they have scheduled an even earlier event, $`t_{\mathrm{step}}`$ will be shortened to that event. The communication is done in the following way: Each process sends $`c^{}`$ and $`t_{\mathrm{cell}}(c^{})`$ as a query to the neighboring processes, which check the adjacent cells for their event times and reply with the minimum thereof. If this answer is smaller, then it is used as $`t_{\mathrm{step}}`$ instead (see Fig. 5, step 1). Periodic tasks like data output, load-balancing or state saving, which require synchronization of the processes can shorten $`t_{\mathrm{step}}`$ even more (see Fig. 5, step 2).

One could think of more rigid policies to determine the length of a step, see e. g. lubachevsky92 . We have also tried other strategies that reduce the number of rollbacks. But on the other hand, they reduce parallelism, too. A large part of the processes are idling, the communication overhead increases, and the overall performance goes down.

After the communication phase, parallel event processing starts and proceeds until the calculated point of time $`t_{\mathrm{step}}`$ (see Fig. 5, step 3), i. e. on average $`𝒪\left(C/C_{\mathrm{border}}\right)`$ iterations are processed. If the algorithm encounters an earlier border zone event which has not been anticipated, the event processing step is stopped immediately to prevent the occurence of a causality error. But normally, the last processed event is a regular border zone event. Then, after the last particle state update, the state information has to be communicated to the neighboring processes.

When a process has finished its computing step, it reads the particle state messages that it has received during this step and adapts its data structures accordingly (see Fig. 5, step 4). Real particles can only be affected as collision partners of virtual particles. For a virtual particle there are several possibilities: A virtual particle can become a real particle by changing from one process to another, it can remain a virtual particle, but with a different position, velocity or cell number, or finally it can emerge on or disappear from the virtual particle list, because it enters or leaves the border zone, respectively. If a virtual particle becomes real and newly calculated events for this particle violate causality, i. e. they are happening before the local simulation time, an error is signaled. Moreover, the processes check the border zone cell times and compare them with the replies to their neighboring processes. If the actual cell time is smaller than the reply (which corresponds to the current time of the neighbor process), a causality error has already happened with high probability, because a collision has been missed on the neighbor process, and an error is signaled, too (see Fig. 5, step 5).

If no error is detected, the simulation continues with the next step. Otherwise a global rollback is performed (see Fig. 5, step 6) and the erroneous simulation step is restarted from the latest saved state (see section III.3). To prevent the reappearance of the same error, a synchronization barrier is issued at the time when the error occured.

Other synchronization barriers can stem from periodic tasks such as state saving, load-balancing or data output (see Fig. 5, step 7).

However, most of all parallel loops do not include synchronization of the processes, but only communication between them.

### III.5 Parallel Performance

First of all, parallelization imposes several performance penalties on the algorithm, among them communication overhead, idle time, state saving, and rollbacks after erroneous computation steps. On the other hand, the computational costs are shared among several processes. If the latter outweighs the former, parallelization can be considered succesful. The same holds true for memory requirements. State saving makes a copy of the whole system state and therefore doubles the memory usage. But here as well parallelization combines the limited resources of single processes.

Figs. 7-7 show that the speed-up of the parallelization (i. e. the reciprocal of the simulation time) for fixed system size is approximately proportional to $`P^{1/2}`$, where $`P`$ is the number of processes.

Furthermore, as shown in Fig. 8, the parallelization is also scalable if the number of processes is chosen proportional to the system size. There are deviations for small $`P`$ when the effect of the parallelization overhead is large. In addition, it has been shown in lubachevsky92  that the error recovery method presented in section III.3 is not scalable for $`P\mathrm{}`$. But for practical purposes, at least until $`P=128`$, the algorithm remains scalable.

We have tested the parallel algorithm on a computer cluster with Pentium III 650 MHz processors and a 1.28 GBit/s switched LAN. On a supercomputer with optimized communication hardware, scalability should be even better.

For real physical applications with more than $`10^2`$ collisions per particle, the maximal number of particles typically is limited by both available memory (see Fig. 9) and computing time to about $`10^6`$ particles per process, i. e. a total number of particles of about $`10^8`$.

## IV Summary

We have demonstrated how to parallelize event-driven molecular dynamics successfully. Our algorithm is based on some ideas from lubachevsky92 , i. e. we have used an optimistic parallelization approach which performs a rollback protocol if a causality error occurs. But the algorithm is enhanced in several ways:

Firstly, we have implemented dynamic load-balancing, which makes simulation of inhomogeneous systems possible. Computing time is further reduced by an adaptive linked-cell structure which determines the optimal cell sizes.

Secondly, we have transferred the shared memory approach of lubachevsky92  to distributed memory architectures as well. The parallelization has been realized with MPI. In order to minimize idle time, we have made use of asynchronous communication, i. e. send and receive actions are decoupled from each other. In addition, the amount of communication is limited to a minimum: The event processing can continue steadily until a border zone event takes place on the local process or is expected to take place on a neighboring process. Only then the calculation has to be interrupted in order to communicate with the neigboring processes. So, even on a cluster with rather poor communication hardware, parallelization yields a speed-up proportional to $`P^{1/2}`$ — at least up to $`P=128`$ parallel processes.

Thirdly, we have optimized the data structure of the algorithm. With event-driven molecular dynamics insufficient memory is often a more serious problem than computing time. In total our optimizations reduce memory requirements to one third as compared to the method proposed in lubachevsky92 . This enabled us to perform simulations of real physical problems with up to $`10^8`$ particles.


text:
# Theoretical uncertainties in sparticle mass predictions from computational tools

## 1 Introduction

Weak-scale supersymmetry (SUSY) is well motivated  because it solves the technical hierarchy problem, removing ultra-violet quadratic divergent corrections to the Higgs mass. Viable weak-scale supersymmetry implies over a hundred “soft breaking” terms in the MSSM Lagrangian at the weak scale, making general analysis intractable. By making model-dependent assumptions on the origin of SUSY breaking for the MSSM fields, one often provides relations between the weak-scale SUSY breaking parameters, vastly reducing the available parameter space. Many theoretical schemes of SUSY breaking exist, the most popular of which can be classified by the mechanism that mediates SUSY breaking from a hidden sector to the MSSM fields. In this paper, we will use minimal supergravity (mSUGRA) , minimal anomaly mediation (mAMSB)  and minimal gauge (mGMSB)  mediation.

Supposing SUSY is discovered at a present or future collider, it will be a major challenge to measure the SUSY breaking parameters with good accuracy. In order to determine the free parameters, physical sparticle masses (or kinematical variables related to them) will be measured. These sparticle masses must be turned into SUSY breaking parameters, typically at some high scale. Provided enough information is collected, this will allow tests on the relations between the high-scale parameters. For example, a simple question already addressed  is: do the gaugino masses unify, and if so, do they unify at the same scale as the gauge couplings? A double affirmative would be strong support in favour of SUSY GUT-type schemes. On the other hand, there are many other possibilities: for example, string models can be non-universal , or unify at an intermediate scale.

There is a large literature upon expected empirical errors on the observables both a Linear Collider  and the LHC . It will be essential to know the theoretical uncertainty if we are to discriminate models of SUSY breaking. There are several studies of how empirical errors propagate into uncertainties on fundamental-scale SUSY breaking parameters in the literature .

The deduction of high-scale SUSY breaking parameters from observables inevitably involves theoretical errors coming from the level of approximation used (e.g. neglected higher order terms), and it is these uncertainties that we study here. We use four modern codes to calculate MSSM spectra: ISAJET 7.64 , SOFTSUSY 1.71 , SPHENO 2.0  and SUSPECT 2.101 . We use the differences in results between the codes to define the current theoretical uncertainties. We have done our best to eliminate differences due to bugs by examining the relevant parts of codes in detail if there was an obvious large discrepancy. However, it would be unrealistic to claim that all of the codes are completely bug free. We therefore take a practical interpretation of ‘theoretical uncertainty’: after all, when fitting experimental data to SUSY breaking models, one must use one of the available computational tools. The precise implementation of the known higher-order corrections differs and has been found to produce significantly different results (for example, using different scales for parameters in the highest-order corrections). Therefore, certainly the differences in results between the codes are due (at least to a large part) to unknown higher-order corrections, matching more traditional notions of ‘theoretical uncertainty’.

Previously, there have been some studies of differences between various calculations of the MSSM spectrum. 2-5$`\%`$ differences in sparticle masses were noticed  between various particles along the Snowmass (SPS)  model lines 1a and 6 between the SOFTSUSY 1.2, ISAJET 7.51 and SUSPECT 2.0 programs. Model line 2 (the so-called focus point line, with very heavy scalar sparticles) was observed to show huge 30% differences in the masses of the weak gauginos. Differences in the spectra and branching ratios were also observed in  comparing ISAJET 7.58, SUSYGEN 3.00 and PYTHIA 6.2. Moreover, large 10$`\%`$ level differences between SOFTSUSY 1.3, ISAJET 7.58 and feynSSG were observed  in the mSUGRA Post-LEP benchmarks . Some of the benchmarks in some of the codes were found to not be consistent with electroweak symmetry breaking, unless one fiddled with the Standard Model inputs $`m_t`$ and $`\alpha _s(M_Z)`$. These were noticeably points E,F (focus-points) and K,M (high $`\mathrm{tan}\beta `$ and high $`m_0,m_{1/2}`$ points). Likewise, the input parameters had to be adjusted in  to get similar spectra from ISAJET 7.51 and SSARD. Initial results highlighting the differences between the predictions in the focus-point and high $`\mathrm{tan}\beta `$ regimes have already been presented by the current authors as a conference proceeding . We will present the main results of this last work here for completeness, updated to state-of-the art calculations.

In this paper, we push these initial observations of theoretical uncertainties further by (i) comparing them to the expected experimental accuracies at the LHC and a future LC and (ii) examining their effect upon future empirical fits to fundamental-scale SUSY breaking schemes. We take expected statistical errors upon sparticle observables from previous studies of mSUGRA at the LHC  at two of the LHC benchmark points. We then perform a fit to mSUGRA at each point with each of the four codes. The statistical precision of the fitted fundamental mSUGRA parameters may then be compared with the theoretical error by looking at the differences between the four fits. We next quantify and present the theoretical error in the coloured sparticle masses at the SPS points. We also re-examine the SPS points providing theoretical uncertainties on mass predictions for sparticles that may kinematically be accessed at a 500 GeV LC. Quantification of errors at particular points will give us an idea of their magnitude and whether or not they significantly depend upon the SUSY breaking scheme (or point).

In section 2, we introduce the four state-of-the art calculations we will use and their level of approximation: ISAJET 7.64, SOFTSUSY 1.71, SPHENO 2.0 and SUSPECT 2.101. ‘Tricky’ regions of the MSSM parameter space, where it is difficult to make accurate predictions, are discussed in section 3: large $`\mathrm{tan}\beta `$ and the focus-point regime. In section 4.1, we perform the LHC mSUGRA empirical fits and then quantify the theoretical error on squark and gluino masses for the SPS points. In section 4.2, we quantify theoretical errors on masses relevant for a 500 GeV Linear Collider. Finally, there are conclusions and an outlook in section 5.

## 2 The codes

We compare the latest versions of four public SUSY renormalisation group evolution (RGE) codes which we think constitute a representative sample of such programs: ISAJET 7.64, SOFTSUSY 1.71, SPHENO 2.0 and SUSPECT 2.101. The basic principle of the SUSY mass spectrum calculation is the same in all programs: Gauge and Yukawa couplings are taken as input parameters at the electroweak scale. However, in the MSSM, in order to define them in the $`\overline{\mathrm{DR}}`$ scheme, from experimental observables, one must first subtract threshold corrections from sparticles. The sparticle spectrum is unknown at this stage and so to begin the calculation, some guess is made for the soft SUSY breaking parameters and spectrum. The MSSM parameters are then run to the high scale $`M_X`$ by RGEs. At $`M_X`$, boundary conditions are imposed on the SUSY breaking parameters. Couplings and SUSY parameters are then run back down to $`M_{SUSY}\sqrt{m_{\stackrel{~}{t}_1}m_{\stackrel{~}{t}_2}}`$ where the $`\mu `$ and $`B`$ MSSM Higgs potential parameters are set in order to give correct radiative electroweak symmetry breaking consistent with an input value of $`\mathrm{tan}\beta =v_2/v_1`$ ($`v_1`$, $`v_2`$ being the two Higgs fields’ vacuum expectation values (VEVs)). The SUSY masses are calculated and radiative corrections are applied, and the parameters are run down to the electroweak scale. Finally, the whole process is iterated in order to obtain a stable solution.

An overview of which corrections are implemented in each of the four programs is given in Table 1. From this table we already expect some differences in the results due to the different levels of radiative corrections applied. In particular ISAJET 7.64 has no finite radiative corrections to most of the sparticle masses. However, we note that even if each column of the table were identical, one could expect different numerical results from each of the codes. In practise, if a quantity is calculated at one-loop, one has the freedom (at one-loop accuracy) of using whatever scale one desires for parameters in the one-loop correction itself. The difference between using, for instance, pole or running masses or couplings derived in the $`\overline{\mathrm{MS}}`$ or $`\overline{\mathrm{DR}}`$ scheme in one-loop corrections is formally of higher-loop order, but leads to non-negligibly different results. Indeed, one can roughly estimate the effects of higher-loop terms by the difference in predicted masses by varying the scale of parameters in the highest loop included.

Let us now point out some of the differences in more detail. A subtle but important issue is the treatment of Yukawa couplings. In general, they are derived from the quark masses as

$$h_t=\sqrt{2}m_t/v_2,h_b=\sqrt{2}m_b/v_1,$$
(1)
where $`m_t`$ and $`m_b`$ are the running top and bottom quark masses in the $`\overline{\mathrm{DR}}`$ scheme. Obviously, it makes a difference whether the VEVs $`v_{1,2}`$ are running (SOFTSUSY 1.71, SPHENO 2.0, SUSPECT 2.101) or not (ISAJET 7.64). It makes an important difference how the $`\overline{\mathrm{DR}}`$ masses are calculated from the pole or $`\overline{\mathrm{MS}}`$ masses. We first discuss Standard Model threshold corrections, and afterwards the sparticle loop corrections.

$``$ ISAJET 7.64 takes 2-loop QCD corrections to $`m_t`$ into account, including the shift from the $`\overline{\mathrm{MS}}`$ to the $`\overline{\mathrm{DR}}`$ scheme as

$$m_t(M_t)_{\mathrm{SM}}^{\overline{\mathrm{DR}}}=M_t\left[1+\frac{5\alpha _s}{3\pi }+\left(16.111.04\left(5\frac{6.63}{M_t}\right)\right)\left(\frac{\alpha _s}{\pi }\right)^2\right]^1,$$
(2)
where $`M_t`$ is the top pole mass and $`\alpha _s=\alpha _s(M_t)`$ at 3-loops in the $`\overline{\mathrm{MS}}`$ scheme. For the bottom quark mass, ISAJET 7.64 takes a hard-coded value of $`m_b(M_Z)_{\mathrm{SM}}^{\overline{\mathrm{DR}}}=2.82`$ GeV .

$``$ SOFTSUSY 1.71 and SPHENO 2.0 calculate both $`h_t`$ and $`h_b`$ at $`Q=M_Z`$. The $`\overline{\mathrm{DR}}`$ top mass is related to the pole mass by 2-loop QCD :

$$m_t(Q)_{\mathrm{SM}}^{\overline{\mathrm{DR}}}=M_t\left[1\frac{\alpha _s}{3\pi }\left(53L\right)\alpha _s^2\left(0.538\frac{43}{24\pi ^2}L+\frac{3}{8\pi ^2}L^2\right)\right]$$
(3)
where111After the publication of this paper, it was noted that eq. (3) incorrectly describes the 2-loop QCD corrections to the top mass: eq. (3) holds for $`L=\mathrm{ln}(m_t^2(Q)/Q^2)`$, not for $`L=\mathrm{ln}(M_t^2/Q^2)`$. SOFTSUSY 1.71 and SPHENO 2.0 therefore contain the incorrect formula, which ought to be
$$m_t(Q)_{\mathrm{SM}}^{\overline{\mathrm{DR}}}=M_t\left[1\frac{\alpha _s}{3\pi }\left(53L\right)\alpha _s^2\left(0.876\frac{91}{24\pi ^2}L+\frac{3}{8\pi ^2}L^2\right)\right]$$
(2.3a)  for $`L=\mathrm{ln}(M_t^2/Q^2)`$. We thank D. R. T. Jones for drawing our attention to this issue. $`L=\mathrm{ln}(M_t^2/Q^2)`$. For the $`b`$ quark, 3-loop relations  are used to calculate $`m_b(M_b)_{\mathrm{SM}}^{\overline{\mathrm{MS}}}`$ from the $`b`$ pole mass, which is run to $`M_Z`$ by 3-loop RGEs . $`m_b(M_Z)_{\mathrm{SM}}^{\overline{\mathrm{MS}}}`$ is then related to the $`\overline{\mathrm{DR}}`$ mass by

$$m_b(M_Z)_{\mathrm{SM}}^{\overline{\mathrm{DR}}}=m_b(M_Z)_{\mathrm{SM}}^{\overline{\mathrm{MS}}}\left[1\frac{\alpha _s}{3\pi }\frac{35\alpha _s^2}{72\pi ^2}+\frac{3g_2^2}{128\pi ^2}+\frac{13g_1^2}{1152\pi ^2}\right]$$
(4)
in SOFTSUSY 1.71. SPHENO 2.0 neglects the last two terms of eq. (4).

$``$ SUSPECT 2.101 uses 2-loop relations and 2-loop RGEs to derive $`m_b(M_b)_{\mathrm{SM}}^{\overline{\mathrm{MS}}}`$ and $`m_t(M_t)_{\mathrm{SM}}^{\overline{\mathrm{MS}}}`$ from the $`b`$ and $`t`$ pole masses. For the conversion to the $`\overline{\mathrm{DR}}`$ scheme, eqs. (3) and (4) are applied (with $`L=0`$ in eq. (3) since $`Q=M_t`$).
Using $`Q=M_t`$ to define $`m_t(Q)`$ (ISAJET 7.64, SUSPECT 2.101) is in principle more accurate than using $`Q=M_Z`$ (SOFTSUSY 1.71, SPHENO 2.0) since then the QCD logs between $`M_Z`$ and $`M_t`$ are re-summed. We refer the reader to the respective manuals  for more details.

The next step is to include SUSY loop corrections. For $`m_t`$, all four programs apply full 1-loop SUSY corrections according to . For $`m_b`$, the full 1-loop (ISAJET 7.64, SPHENO 2.0) or the leading (SOFTSUSY 1.71, SUSPECT 2.101) SUSY corrections are included. The $`\mathrm{tan}\beta `$ enhanced corrections to $`m_b`$ are re-summed as given in  in all four programs, see . Still, there are some important differences in the implementation of these corrections. For example, the $`\stackrel{~}{g}\stackrel{~}{t}`$ correction to $`m_t`$ is

$`\left({\displaystyle \frac{\mathrm{\Delta }m_t}{m_t}}\right)^{\stackrel{~}{g}\stackrel{~}{t}}`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{3\pi }}\{B_1(m_t^2,m_{\stackrel{~}{g}}^2,m_{\stackrel{~}{t}_1}^2)+B_1(m_t^2,m_{\stackrel{~}{g}}^2,m_{\stackrel{~}{t}_2}^2)`$ (5)
$`+{\displaystyle \frac{m_{\stackrel{~}{g}}}{m_t}}\mathrm{sin}2\theta _{\stackrel{~}{t}}[B_0(m_t^2,m_{\stackrel{~}{g}}^2,m_{\stackrel{~}{t}_1}^2)B_0(m_t^2,m_{\stackrel{~}{g}}^2,m_{\stackrel{~}{t}_2}^2)]\}.`$
ISAJET 7.64 adds this correction at $`m_t`$, but using a scale $`Q=\sqrt{m_{\stackrel{~}{t}_L}m_{\stackrel{~}{t}_R}}`$, (this is done to avoid double counting of logarithmic corrections which are included via step functions in the RGEs) and $`\alpha _s=\alpha _s(m_{\stackrel{~}{g}})`$ in eq. (5), while SOFTSUSY 1.71 and SPHENO 2.0 calculate it at $`Q=M_Z`$ and SUSPECT 2.101 at $`Q=M_t`$. Accordingly, $`m_t=m_t(Q)_{\mathrm{SM}}^{\overline{\mathrm{DR}}}`$ in the term $`m_{\stackrel{~}{g}}/m_t`$ enters with different values in all four programs. Due to differences in the inclusion of finite radiative corrections to sparticle masses, the gluino and stop masses in (5) vary from program to program. In particular, ISAJET 7.64 calculates the corrections to $`m_{\stackrel{~}{g}}`$ at $`Q=m_{\stackrel{~}{g}}`$ while SUSPECT 2.101 calculates them at $`Q=M_Z`$, which leads to quite different gluino masses. Analogous differences exist in the other contributions to $`(\mathrm{\Delta }m_t)^{SUSY}`$ as well as in the calculation of $`(\mathrm{\Delta }m_b)^{SUSY}`$.

Another comment is in order concerning $`\alpha _s`$. The value of $`\alpha _s(M_Z)`$ from experiment is given in the $`\overline{\mathrm{MS}}`$ scheme. SOFTSUSY 1.71, SPHENO 2.0 and SUSPECT 2.101 take $`\alpha _s(M_Z)^{\overline{\mathrm{MS}}}`$ as input and convert it to the $`\overline{\mathrm{DR}}`$ scheme

$$\alpha _s^{\overline{\mathrm{DR}}}(M_Z)=\frac{\alpha _s(M_Z)^{\overline{\mathrm{MS}}}}{1\mathrm{\Delta }\alpha _s},$$
(6)
where

$$\mathrm{\Delta }\alpha _s=\frac{\alpha _s(M_Z)}{2\pi }\left[\frac{1}{2}\frac{2}{3}\mathrm{ln}\left(\frac{M_t}{M_Z}\right)2\mathrm{ln}\left(\frac{m_{\stackrel{~}{g}}}{M_Z}\right)\frac{1}{6}\underset{\stackrel{~}{q}}{}\underset{i=1,2}{}\mathrm{ln}\left(\frac{m_{\stackrel{~}{q}_i}}{M_Z}\right)\right].$$
(7)
In SUSPECT 2.101, the log terms are not added explicitly but included via threshold functions in the RGEs. Also ISAJET 7.64 re-sums the logs in (7) via step-function decoupling in the RGEs. The finite term, however, is not taken into account. The difference due to the finite term is small but relevant (1%). One could in principle interpret the input $`\alpha _s(M_Z)`$ in ISAJET 7.64 as already being the effective Standard Model $`\overline{\mathrm{DR}}`$ value. Since, however, for some corrections ISAJET 7.64 effectively takes hard-wired values of $`\alpha _s(M_Z)^{\overline{\mathrm{MS}}}=0.118`$, we take the canonical input value of $`\alpha _s(M_Z)^{\overline{\mathrm{MS}}}=0.118`$ for all codes.

## 3 Tricky corners of SUSY parameter space

As a general point, when quantifying errors on predicted observables, we assume that the results from the codes follow a Gaussian probability distribution. This is a priori unjustified, but we prefer it to quoting minimum and maximum values because we believe that the additional information included in the variance is desirable. For example, if three codes all provided identical results and one gave a result further away, we think that the true uncertainty ought to be less than the range of minimum-maximum results.

### 3.1 Large $`\mathrm{tan}\beta `$

Large $`\mathrm{tan}\beta `$ has always been recognised as a difficult case since it requires a thorough treatment of the bottom Yukawa coupling. Figure 1 shows $`h_b`$ of the four different programs as a function of $`\mathrm{tan}\beta `$ in the mSUGRA model. We see that SOFTSUSY 1.71 and SPHENO 2.0 agree very well on $`h_b`$, and there is also good agreement with ISAJET 7.64. Comparing only these programs we would assign a $`\stackrel{<}{}\mathrm{\hspace{0.33em}3}\%`$ uncertainty on $`h_b`$ even for very large $`\mathrm{tan}\beta `$. The agreement with SUSPECT 2.101 is however not so good, and we find 4 – 8% uncertainty taking all four programs into account. The effect of re-summing the $`\mathrm{tan}\beta `$ enhanced SUSY loop corrections can be seen when comparing the solid and dotted lines in fig. 1, the solid line being the result of ISAJET 7.64, where the $`(\mathrm{\Delta }m_b)^{SUSY}`$ corrections are resummed, and the dotted one being the result of ISASUSY 7.58, where this re-summation is not applied.

The bottom Yukawa coupling has its largest effect in the Higgs sector when it is large (at high $`\mathrm{tan}\beta `$): the evolution of $`m_{H_1}^2`$ is driven by $`h_b`$,

$$\frac{dm_{H_1}^2}{dt}\frac{3}{8\pi ^2}h_bX_b+\mathrm{},X_b=(m_{\stackrel{~}{Q}}^2+m_{\stackrel{~}{D}}^2+m_{H_1}^2+A_b^2),$$
(8)
where $`t=\mathrm{ln}Q`$, $`Q`$ being the renormalisation scale. Differences in $`m_{H_1}^2`$ directly translate into the physical Higgs boson masses since

$$m_A^2=\frac{1}{c_{2\beta }}\left(\overline{m}_{H_2}^2\overline{m}_{H_1}^2\right)+\frac{s_\beta ^2t_1}{v_1}+\frac{c_\beta ^2t_2}{v_2}M_Z^2.$$
(9)
Here $`\overline{m}_{H_i}^2=m_{H_i}^2t_i/v_i`$, $`i=1,2`$, and $`t_{1,2}`$ are the tadpole contributions. The self energies of $`Z`$ and $`A`$ have been neglected in eq. (9).

Figure 2 shows the Higgs boson masses obtained by the four programs as a function of $`\mathrm{tan}\beta `$ (the results obtained by ISAJET 7.58 are again shown as dotted lines in fig. 2). Let us first discuss the masses of $`A^0`$ and $`H^\pm `$. For $`\mathrm{tan}\beta =10`$ – 50, we find differences in $`m_A`$ and $`m_{H^\pm }`$ of about 10 – 50 GeV, dominated over most of the $`\mathrm{tan}\beta `$ range by the difference between SPHENO 2.0 and SUSPECT 2.101. This has to be compared with differences of 100 GeV and more encountered with earlier versions, see for instance the dotted lines representing ISAJET 7.58. Assuming the error to be Gaussian, we now have $`\mathrm{\Delta }m_{A,H^\pm }\pm 10`$ GeV at $`\mathrm{tan}\beta =25`$ and $`\mathrm{\Delta }m_{A,H^\pm }\pm 20`$ GeV at $`\mathrm{tan}\beta =50`$. The bottom Yukawa coupling is, however, not the only source of differences in $`m_A`$. Another source is, for example, whether one uses running or pole values for masses in the calculation of one-loop tadpoles $`t_{1,2}`$. Also, the scale and scheme of parameters in the one-loop expressions for the tadpoles all vary. These differences are formally of higher order and indeed each program has a different approach.

The situation is somewhat different for the $`(h^0,H^0)`$ system because here additional radiative corrections are necessary. It is well known that these involve a theoretical uncertainty on $`m_{h^0}`$ of about 3 GeV , evidence of which can be seen in fig. 2. For completeness we note that SUSPECT 2.101 offers various choices of Higgs mass calculations. In fig. 2, we have used its default $`m_{h^0}`$ routine, i.e. ichoice(10) = 0, shown as a black dash-dot-dotted line. If we use instead SUSPECT 2.101 with FeynHiggsFast, ichoice(10) = 3, we get $`m_{h^0}115`$ GeV, shown as a grey dash-dot-dotted line in fig. 2. This will be relevant later in this paper when we discuss mSUGRA fits to LHC data.

### 3.2 Focus point

For large $`m_0`$, the running of $`m_{H_2}^2`$ becomes very steep and very sensitive to the top Yukawa coupling:

$$\frac{dm_{H_2}^2}{dt}\frac{3}{8\pi ^2}h_tX_t+\mathrm{},X_t=(m_{\stackrel{~}{Q}}^2+m_{\stackrel{~}{U}}^2+m_{H_2}^2+A_t^2).$$
(10)
As a result, the $`\mu `$ parameter given by

$$\mu ^2=\frac{\overline{m}_{H_1}\overline{m}_{H_2}^2\mathrm{tan}^2\beta }{\mathrm{tan}^2\beta 1}\frac{1}{2}M_Z^2$$
(11)
becomes extremely sensitive to $`h_t`$. This is visualised in fig. 3 where we show $`h_t`$ and $`\mu `$ as functions of $`m_0`$. The other parameters are $`m_{1/2}=300`$ GeV, $`A_0=0`$, $`\mathrm{tan}\beta =10`$ and $`\mu >0`$. There is reasonable agreement on $`\mu `$ up to $`m_01`$ TeV, and the differences observed for $`m_0\stackrel{<}{}\mathrm{\hspace{0.33em}2}`$ TeV are phenomenologically not so important. For larger values of $`m_0`$ we observe, however, large discrepancies between the four programs. These lead to completely different chargino/neutralino properties for very large $`m_0`$, and likewise to very different excluded regions, depending on which program is used. The situation is, however, already much better than the one reported in , c.f. the dashed lines of ISASUSY 7.58 in fig. 3 (earlier versions of the other codes also gave quite different results).

In order to understand the behaviour of $`\mu `$ in fig. 3 it is useful to write eq. (11) in the form

$$\mu ^2c_1m_0^2+c_2m_{1/2}^20.5M_Z^2.$$
(12)
Approximate analytical expressions for $`c_1`$ and $`c_2`$ can be found e.g., in . For $`A_0=0`$ and $`\mathrm{tan}\beta =10`$ we get

$$c_1\left(\frac{m_t}{156.5\mathrm{GeV}}\right)^21,c_2\left(\frac{m_t}{102.5\mathrm{GeV}}\right)^20.52.$$
(13)
Figure 4a shows contours of constant $`\mu `$ in this approximation in the $`(m_0,m_t)`$ plane. Notice the fast increasing dependence on $`m_t`$ for increasing $`m_0`$. For $`m_t156`$ – 157 GeV, $`\mu `$ becomes almost independent of $`m_0`$. This is a signal of the actual focus point behaviour , which is defined as the insensitivity of $`m_{H_2}`$ to its GUT scale value. Eq. 11 then implies that $`\mu `$ is insensitive to $`m_0`$, since $`m_{H_1}`$ appears with a suppression factor of $`\mathrm{tan}^2\beta 1`$. Interpreting $`m_t`$ in eq. (13) as $`m_t(M_{SUSY})=h_t(M_{SUSY})v_2(M_{SUSY})/\sqrt{2}`$, we can directly relate the $`m_0`$ dependence of $`\mu `$ in fig. 3 to that of $`h_t`$.

Some more comments are in order. Firstly, SUSPECT 2.101 has only 1–loop RGEs for scalar SUSY parameters. For large $`m_0`$, the 2–loop terms lead to $`𝒪(10\%)`$ correction and should thus be taken into account. Secondly, as already mentioned in section 2, the sparticle masses that enter the radiative corrections have different values in different codes. In particular, there are large differences in the gluino masses between SUSPECT 2.101 and the other codes, as illustrated in fig. 4b. Since the gluino mass enters $`h_t`$ via the $`\stackrel{~}{g}\stackrel{~}{t}`$ correction eq. (5), this may account for the different slope of $`h_t`$ (and consequently of $`\mu `$) in SUSPECT 2.101 compared to the other programs, as evident in fig. 3.

It thus seems that for reliable results in the large $`m_0`$ region, a more complete calculation of the top Yukawa coupling at the 2–loop level is necessary. This should include finite radiative corrections to all sparticle masses at the full 1–loop level.

## 4 Comparison of theoretical and experimental uncertainties

### 4.1 Fits of mSUGRA parameters to LHC data

In the ATLAS TDR , a case study was made of fitting the mSUGRA model to possible measurements of six reference scenarios. We have re-analysed these fits for two of these scenarios, LHC Point 1 and Point 2. Here squarks and gluinos are produced with the dominant decays $`\stackrel{~}{g}q\stackrel{~}{q}_{L,R}`$, $`\stackrel{~}{q}_L\stackrel{~}{\chi }_2^0q\stackrel{~}{\chi }_1^0hq`$, $`\stackrel{~}{q}_R\stackrel{~}{\chi }_1^0q`$. The assumed measurements for the two points for low ($`=30`$ fb-1) and high ($`=300`$ fb-1) luminosity are given in table 2. They were estimated in ref.  by using ISAJET 7.34 and simulating the ATLAS experiment to determine expected empirical errors.222It is beyond the scope of the present paper to re-perform the experimental analysis in order to have the numbers in table 2 more up-to-date.

With each of the programs under discussion we have performed a $`\chi ^2`$ fit of the mSUGRA parameters $`m_0`$, $`m_{1/2}`$ and $`\mathrm{tan}\beta `$ to the data of table 2, taking $`A_0=0`$ and $`\mu >0`$. The results are listed in tables 3 and 4. The quoted errors are at 1$`\sigma `$ (68.3% C.L.) from a simultaneous fit of all three parameters, i.e. $`\mathrm{\Delta }\chi ^2=3.53`$. In case of SUSPECT 2.101, we have used its default option, ichoice(10) = 0, for the calculation of the $`h^0`$ mass. When linking SUSPECT 2.101 with FeynHiggsFast, ichoice(10) = 3, the results for $`m_0`$ and $`m_{1/2}`$ practically do not change. However, we get much lower values for $`\mathrm{tan}\beta `$: $`\mathrm{tan}\beta =1.7\pm 0.1`$ for Point 1, and $`\mathrm{tan}\beta =6.06\pm 2.06`$ for Point 2 at high luminosity.

It is interesting to note that not only the central values but also the size of the errors can be quite different. Similarly, the minimum $`\chi ^2`$, which is a measure of the quality of the fit, can show large variations. To make the comparison easier and to visualise correlations and non-Gaussian effects, we show in Figs. 5 and 6 contours of 68% and 95% C.L. for the fits of Points 1 and 2 in the $`(m_0,m_{1/2})`$ and $`(m_0,\mathrm{tan}\beta )`$ planes. The values of $`\mathrm{\Delta }\chi ^2`$ used for these confidence levels are based upon a simultaneous two-parameter fit, i.e. $`\mathrm{\Delta }\chi ^2=2.3`$ and 5.99. The third parameter, $`\mathrm{tan}\beta `$ or $`m_{1/2}`$, is always fixed to its best-fit value. As one can see, the error ellipses have only little or even no overlap. We therefore conclude that the theoretical uncertainty is about the same size as the statistical one in the fitted quantities.

One might expect that the main source of these differences is the theoretical uncertainty on $`m_h`$. This is indeed the case for the determination of $`\mathrm{tan}\beta `$. However, $`m_h`$ has only little influence on the fit of $`m_0`$ and $`m_{1/2}`$. This becomes clear when using e.g., SUSPECT 2.101 with different routines for the Higgs mass calculation. In fig. 6a, we show the results of SUSPECT 2.101 \+ FeynHiggsFast, ichoice(10) = 3, as dashed contours in addition to those obtained with its default Higgs mass routine (ichoice(10) = 0, solid contours). In fig. 6b, we have omitted the default SUSPECT 2.101 results to avoid confusion of the many lines. They would look similar to the SPHENO 2.0 contours but centred at $`m_0=450`$ GeV and without an upper limit on $`\mathrm{tan}\beta `$. We thus conclude that the uncertainties in $`m_0`$ and $`m_{1/2}`$ mainly come from the differences in the programs pointed out in Sects. 23, and not from different Higgs mass calculations.

We next address the question of how the theoretical uncertainty depends on the SUSY parameter point. Copious quantities of squarks and gluinos are expected to be produced at the LHC, leading to a fairly precise measurement of their masses, particularly if $`m_{\chi _1^0}`$ is determined accurately by a LC. In table 5 we compare the $`\stackrel{~}{g}`$, $`\stackrel{~}{u}_L`$, $`\stackrel{~}{u}_R`$ and $`\stackrel{~}{t}_1`$ masses obtained by the four programs for the Snowmass (SPS) points . Assuming a Gaussian distribution, we quote the variance of these masses as the theoretical error, i.e. $`\delta m_X=\sqrt{\frac{1}{N1}_i[(m_X)_i\overline{m}_X]^2}`$ where $`\overline{m}_X`$ is the mean of $`(m_X)_i`$ and $`N=4`$ in our case. We make the following observations: (i) the absolute theoretical uncertainty (in GeV) varies from point to point; (ii) the typical relative uncertainty in mSUGRA and mGMSB scenarios in generic (i.e. not tricky) regions of parameter space is about 2 – 5%; (iii) in some cases, in particular in focus point and mAMSB scenarios, the relative uncertainty is larger, about 5 – 10%; (iv) in any case, the theoretical error is of the same order of magnitude as the experimental one.

### 4.2 Linear Collider measurements

At a high-luminosity $`e^+e^{}`$ Linear Collider, one expects to measure chargino, neutralino and slepton masses with accuracies at the per-cent or even per-mill level . We thus take the differences in these masses accessible at $`\sqrt{s}=500\mathrm{GeV}`$ as a measure of the theoretical uncertainty. Tables 6 and 7 compare masses obtained by the four programs for the various SPS points. Only those masses kinematically accessible at a 500 GeV LC are shown. The error is again defined as the variance of the masses as in the previous section.

It turns out that the uncertainty in the LSP mass in mSUGRA and mGMSB scenarios is typically a few hundred MeV, depending on the parameter point. An exception is the focus point scenario (SPS2) where $`\delta m_{\stackrel{~}{\chi }_1^0}=1.4`$ GeV. For $`\stackrel{~}{\chi }_2^0`$ and $`\stackrel{~}{\chi }_1^\pm `$, we find uncertainties of about 1 – 3 GeV in mSUGRA and mGMSB scenarios. (Here note that with earlier program versions, especially with ISASUSY 7.51 – 7.63, we had discrepancies of 50% and more for focus point scenarios.) For the sleptons we find typical uncertainties of 1 – 2 GeV. We note that for SPS5, the lighter stop would be accessible, with $`m_{\stackrel{~}{t}_1}=235.6`$, 232.3, 248.4 and 242.6 GeV for ISAJET 7.64, SOFTSUSY 1.71, SPHENO 2.0 and SUSPECT 2.101, respectively, corresponding to an error of 7 GeV, c.f. table 5. In the mAMSB scenario, SPS 9, we have much larger uncertainties of $``$ 8 GeV for $`m_{\stackrel{~}{\chi }_1^0}`$ and $`m_{\stackrel{~}{\chi }_1^+}`$. This is due to the fact that SOFTSUSY 1.71 has 2-loop GUT-scale boundary conditions for mAMSB, while the other programs have only 1-loop boundary conditions. If we enforce 1-loop boundary conditions for all four programs the error decreases to 3 GeV. We therefore expect that if all programs were to use 2-loop boundary conditions the error would be better than 3 GeV.

By comparing table 5 with tables 6, 7 we see that the present theoretical uncertainty in mass predictions is significantly smaller for weakly interacting sparticles than for strongly interacting ones. This is expected, since the former have smaller threshold corrections. However, the errors in tables 6, 7 are still larger by up to an order of magnitude than the expected experimental accuracies at an $`e^+e^{}`$ Linear Collider. Moreover, the differences in the masses produce cross sections and branching ratios that differ by a few per-cent. While this is not a problem for determining SUSY-breaking parameters at the low scale, it will be relevant when relating them to GUT-scale parameters in order to test their unification and to determine the sources of SUSY breaking.

## 5 Conclusions

If sparticles are detected at future colliders, measurements of their properties (and confirmation that they are in fact sparticles) will take place. The question following from this line of investigation will be: what can we learn about SUSY breaking from these measurements? Is there a unification of certain SUSY breaking parameters, and is it possible to distinguish various SUSY breaking models?

Vital ingredients for answers to these questions will be the amount and precision of empirical measurements made at the particular colliders . Another vital ingredient will be the precision with which we can relate the experimentally observed quantities at the TeV scale with the fundamental physics at the high-energy scale. Accurate extrapolations require multi-loop results for the RGEs and the related threshold corrections. In this paper, we have addressed the question: what is the current theoretical uncertainty associated with determining fundamental-scale SUSY parameters? To this end, we compared four public state-of-the-art MSSM spectrum calculations: ISAJET 7.64, SOFTSUSY 1.71, SPHENO 2.0 and SUSPECT 2.101, taking the spread of their results as a measure of the to-date uncertainty. Although this does not correspond to the usual notion of ‘theoretical uncertainty’, it is pragmatic in the sense that (at least) one of the available calculational tools will be used to perform fits if and when the relevant data arrives. The uncertainty was shown to be largest in certain tricky corners of parameter space: the focus point region and high $`\mathrm{tan}\beta `$. However, even in these regions, comparison with previous versions of the codes  shows that the theoretical uncertainty has significantly improved. Sparticle masses in these regions are particularly sensitive to the values of the Yukawa couplings (especially the top Yukawa for the focus point, and the bottom Yukawa for the high $`\mathrm{tan}\beta `$ regime). Slightly different treatments of top and bottom masses can lead to large differences in mass predictions. It is therefore critical that the accuracy in deriving the running top and bottom $`\overline{\mathrm{DR}}`$ Yukawa couplings is at a maximum in any calculation.

We used previous LHC estimates of expected empirical errors at two benchmark points to perform separate fits to mSUGRA with the four different codes. The parameters resulting from the four calculations show a difference comparable to the statistical error upon them, showing that theoretical uncertainties must be taken into account. We then went on to quantify the theory uncertainty associated with squark and gluino masses for the SPS benchmark points. The theory uncertainty was also quantified for parts of the MSSM spectrum that would be kinematically accessible to a 500 GeV Linear Collider for the SPS points. The tables of spectra listed in the present article may be used for further comparisons with other (possibly new or private) computations. The theoretical uncertainties for LHC and Linear Collider observables were shown to significantly depend upon the point in parameter space being considered. It will not be very precise, therefore, to use the theory uncertainties calculated here when analysing actual data, since the SUSY breaking parameter point is a priori unknown.

If sparticles are measured, the correct modus operandi is clear when analysing data: any fits to a fundamental SUSY breaking model should be performed with several state-of-the-art spectrum calculations in order to deduce the level of theoretical uncertainty, which must then be included in final quoted errors. The theory errors must clearly be included when discussing the accuracy required to distinguish between SUSY models on sparticle observables, such as in ref. .

We note that each of the codes used assumed a pure MSSM desert up until the GUT scale $`10^{16}`$ GeV. We know already  that this assumption is called into question because the gauge couplings themselves require $`𝒪(1\%)`$ corrections at the high scale in order to make them unify properly. It is thought that this level of correction is perfectly reasonable in, for example, SUSY GUTs since model-dependent GUT-scale threshold corrections  are expected from (for example) heavy coloured triplets. If the combined theoretical and empirical accuracy is significantly better than $`1\%`$, then observables at colliders could be used to measure these threshold corrections.

Fortunately, theoretical errors in sparticle mass predictions are certainly not static. There has been much progress reducing them recently (especially in the Higgs and electroweak symmetry breaking sectors, see for example ). We expect this trend to continue, which, as our present results indicate, is desirable if we are to disentangle SUSY breaking from experimental observables.

###### Acknowledgments.

We thank Howie Baer, Abdelhak Djouadi, Fabiola Gianotti, Frederick James, Konstantin Matchev and Giacomo Polesello for useful discussions. This work was partly funded by CNRS. W.P. is supported by the Swiss ‘Nationalfonds’.