text:
# Classification of mixed three–qubit states

## Abstract

We introduce a classification of mixed three–qubit states, in which we define the classes of separable, biseparable, W– and GHZ–states. These classes are successively embedded into each other. We show that contrary to pure W–type states, the mixed $`W`$–class is not of measure zero. We construct witness operators that detect the class of a mixed state. We discuss the conjecture that all entangled states with positive partial transpose (PPTES) belong to the $`W`$–class. Finally, we present a new family of PPTES “edge” states with maximal ranks.

The rapidly increasing interest in quantum information processing has motivated the detailed study of entanglement. Whereas entanglement of pure bipartite systems is well understood, the classification of mixed states according to the degree and character of their entanglement is still a matter of intensive research (see ). It was soon realised, that the entanglement of pure tripartite quantum states is not a trivial extension of the entanglement of bipartite systems. Recently, the first results concerning the entanglement of pure tripartite systems have been achieved. There, the main goal has been to generalize the concept of the Schmidt decomposition to three-party systems , and to distinguish classes of locally inequivalent states. The knowledge of mixed tripartite entanglement is much less advanced (see, however, ).

In this Letter we introduce a classification of the whole space of mixed three–qubit states into different entanglement classes. We provide a method to determine to which class a given state belongs (tripartite witnesses). We also discuss the characterization of entangled states that are positive under partial transposition (PPTES). Finally, we introduce a new family of PPTES for mixed tripartite qubits.

Our proposal to classify mixed tripartite–qubit states is done by specifying compact convex subsets of the space of all states, which are embedded into each other. This idea vaguely resembles the classification of bipartite systems by their Schmidt number . However, as shown later our classification does not follow the Schmidt number. Also in this respect, entanglement of tripartite systems differs genuinely from the one of bipartite quantum systems.

Before presenting our results concerning mixed states, we briefly review some of the recent results on pure three–qubit states. Any three–qubit vector (pure state) can be written as

$`|\psi _{GHZ}`$ $`=`$ $`\lambda _0|\mathrm{\hspace{0.17em}000}+\lambda _1e^{i\theta }|\mathrm{\hspace{0.17em}100}+\lambda _2|\mathrm{\hspace{0.17em}101}`$ (1)
$`+`$ $`\lambda _3|\mathrm{\hspace{0.17em}110}+\lambda _4|\mathrm{\hspace{0.17em}111},`$ (2)
where $`\lambda _i0`$, $`_i\lambda _i^2=1`$, $`\theta [0,\pi ]`$ , and $`\{|0,|1\}`$ denotes an orthonormal basis in Alice’s, Bob’s and Charlie’s space, respectively. Apart from separable and biseparable pure states, there exist also two different types of locally inequivalent entangled vectors; the so-called GHZ–type  and W–type. Vectors belonging to GHZ– and W–types cannot be transformed into each other by local operations and classical communication (LOCC). Generically, a vector described by Eq.(2) is of the GHZ–type, while W–vectors can be written as

$$|\psi _W=\lambda _0|\mathrm{\hspace{0.17em}000}+\lambda _1|\mathrm{\hspace{0.17em}100}+\lambda _2|\mathrm{\hspace{0.17em}101}+\lambda _3|\mathrm{\hspace{0.17em}110}.$$
(3)
W–vectors form a set of measure zero among all pure states. Also, given a W–vector one can always find a GHZ–vector as close to it as desired by adding an infinitesimal $`\lambda _4`$–term to the RHS of Eq.(3) . Furthermore, the so-called tangle, $`\tau `$, introduced in , can be used to detect the type, since $`\tau (|\psi _W)=0`$ .

Mixed states of three–qubit systems can be classified generalizing the classification of pure states. To this aim we define (see Fig.1):
$``$ the class $`S`$ of separable states, i.e. those that can be expressed as a convex sum of projectors onto product vectors;
$``$ the class $`B`$ of biseparable states, i.e. those that can be expressed as a convex sum of projectors onto product and bipartite entangled vectors (A–BC, B–AC and C–AB);
$``$ the class $`W`$ of W–states, i.e. those that can be expressed as a convex sum of projectors onto product, biseparable and W–type vectors;
$``$ the class $`GHZ`$ of GHZ–states, i.e. the set of all physical states.
All these sets are convex and compact, and satisfy $`SBWGHZ`$. States in $`S`$ are not entangled. No genuine three–party entanglement is needed to prepare entangled states in the subset $`BS`$. The formation of entangled states in $`WB`$ requires W–type vectors with three–party entanglement, but zero tangle, which is an entanglement monotone decreasing under LOCC . Finally, the class $`GHZ`$ contains all types of entanglement, and in particular, GHZ–type vectors are needed to prepare states from $`GHZW`$. The introduced classes are invariant under local unitary or invertible non-unitary operations, while local POVM’s  can only transform states from a “higher” to a “lower” class.

Notice that since GHZ–vectors can be expressed as the sum of only two product vectors, i.e. $`|GHZ=(|\mathrm{\hspace{0.17em}000}+|\mathrm{\hspace{0.17em}111})/\sqrt{2}`$, whereas the minimum number of product terms forming a W–vector is three , as in the state $`|W=(|\mathrm{\hspace{0.17em}100}+|\mathrm{\hspace{0.17em}010}+|\mathrm{\hspace{0.17em}001})/\sqrt{3}`$, our scheme may seem somehow counterintuitive. Indeed, for bipartite systems, states with lower Schmidt number, i.e. lower number of product terms in the Schmidt decomposition, are embedded into the set of states with higher Schmidt number . One is tempted to extend this classification to tripartite systems as $`SBGHZW`$, where now $`W`$ is the set of all states. However, such generalization is evidently wrong, because the the set of GHZ–states in such classification cannot be closed .

Having established the structure of the set of mixed three–qubit states, we show how to determine to which class a given state $`\rho `$ belongs. To this aim, we use the approach developed previously in the construction and optimisation of witness operators .

We denote the range of $`\rho `$ by $`R(\rho )`$, its rank by $`r(\rho )`$, its kernel by $`K(\rho )`$, and the dimension of $`K(\rho )`$ by $`k(\rho )`$. Following the approach of the best separable approximation (BSA), one can decompose any state $`\rho `$ as a convex combination of a $`W`$–class state and a remainder $`\delta `$,

$$\rho =\lambda _W\rho _W+(1\lambda _W)\delta ,$$
(4)
where $`0\lambda _W1`$, and $`R(\delta )`$ does not contain any W–vector. Maximization of $`\lambda _W`$ leads to the best W–approximation of $`\rho `$. Notice that only for $`\rho `$ belonging to the $`GHZW`$–class, this decomposition is non-trivial, i.e. $`\lambda _W1`$. Also, $`r(\delta )=1`$, since any subspace spanned by two linearly independent GHZ–vectors contains at least one pure state with zero tangle. In fact, given $`|\psi _1`$ and $`|\psi _2`$ with $`\tau (|\psi _1)`$ and $`\tau (|\psi _2)`$ not equal zero, it is always possible to find some $`\stackrel{~}{\alpha },\stackrel{~}{\beta }`$ such that $`|\psi (\stackrel{~}{\alpha },\stackrel{~}{\beta })=\stackrel{~}{\alpha }|\psi _1+\stackrel{~}{\beta }|\psi _2`$ is normalized, and its tangle is zero. Therefore, any W–approximation must have the form:

$$\rho =\lambda _W\rho _W+(1\lambda _W)|\psi _{GHZ}\psi _{GHZ}|.$$
(5)
Similarly, one can express $`\rho `$ in the best biseparable approximation as:

$$\rho =\lambda _B\rho _B+(1\lambda _B)\delta ,$$
(6)
where now $`R(\delta )`$ must not contain any biseparable states, i.e. $`r(\delta )<4`$, since any $`N`$–dimensional subspace of the $`2\times N`$ space contains at least one product vector .

We use the above decompositions to construct operators that detect the desired subset (see ). In analogy to entanglement witnesses and Schmidt witnesses we term these operators tripartite witnesses. The existence of witness operators is a consequence of the Hahn-Banach theorem, which states that a point outside a convex compact set is separated from that set by a hyper-plane. The equation Tr$`(𝒲\rho )=0`$ describes such a hyper-plane, and one calls $`𝒲`$ a witness operator. For example, in our setting, a W–witness is an operator $`𝒲_W`$ such that Tr$`(𝒲_W\rho _B)0`$ holds $`\rho _BB`$, but for which there exists a $`\rho _WWB`$ such that Tr$`(𝒲_W\rho _W)<0`$.

Any GHZ–witness (W–witness) has the canonical form $`𝒲=Qϵ1𝐥`$, where $`Q`$ is a positive operator which has no W–type (B–type) vectors in its kernel; thus $`k(Q)=1`$ ($`k(Q)<4`$). An example of a GHZ–witness is

$$𝒲_{GHZ}=\frac{3}{4}1𝐥P_{GHZ},$$
(7)
where $`P_{GHZ}`$ is the projector onto $`|GHZ`$. The value 3/4 corresponds to the maximal squared overlap between $`|GHZ`$ and a W–vector. This construction guarantees that Tr$`(𝒲_{GHZ}\rho _W)0`$ for any W–state, and since Tr$`(𝒲_{GHZ}P_{GHZ})<0`$, there is a GHZ $``$ W–state which is detected by $`𝒲_{GHZ}`$. The maximal overlap is obtained as follows: due to the symmetry of $`|GHZ`$ we only need to consider W–vectors that are symmetric under the exchange of any of the three qubits. Therefore, we have to consider all local trilateral rotations of $`|\psi _W=\kappa _0|\mathrm{\hspace{0.17em}000}+\kappa _1(|\mathrm{\hspace{0.17em}100}+|\mathrm{\hspace{0.17em}010}+|\mathrm{\hspace{0.17em}001})`$, where $`\kappa _0,\kappa _1`$ are real and $`\kappa _0^2+3\kappa _1^2=1`$. Due to the symmetry, such rotations can be parametrised for all parties as $`|\mathrm{\hspace{0.17em}0}\alpha |\mathrm{\hspace{0.17em}0}+\beta |\mathrm{\hspace{0.17em}1}`$, $`|\mathrm{\hspace{0.17em}1}\beta ^{}|\mathrm{\hspace{0.17em}0}\alpha ^{}|\mathrm{\hspace{0.17em}1}`$, with $`|\alpha |^2+|\beta |^2=1`$. Thus, the overlap $`GHZ|\psi _W`$ is a function of six parameters with two constraints, and can be maximized using Lagrange multipliers. An optimal choice of parameters is $`\kappa _0=0`$, $`\kappa _1=1/\sqrt{3}`$, and $`\beta =\alpha =1/\sqrt{2}`$. This leads to $`|GHZ|\psi _W|_{max}^2=3/4`$.

Analogously, we can construct a W–witness as

$$𝒲_{W_1}=\frac{2}{3}1𝐥P_W,$$
(8)
where $`P_W`$ is now the projector onto a vector $`|W`$, and 2/3 corresponds to the maximal squared overlap between $`|W`$ and a B–vector. Another example of a W–witness is

$$𝒲_{W_2}=\frac{1}{2}1𝐥P_{GHZ},$$
(9)
where now 1/2 is the maximal squared overlap between $`|GHZ`$ and a B–type vector . The W–vector that has maximal overlap with $`|GHZ`$ is detected by $`𝒲_{W_2}`$.

The tripartite witness $`𝒲_{W_2}`$ allows to prove that the class of mixed W $``$ B–states is not of measure zero: consider the family of states in $`𝒞^2𝒞^2𝒞^2`$ given by the convex sum of the identity and a projector onto a W–state,

$$\rho =\frac{1p}{8}1𝐥+pP_W.$$
(10)
Obviously, the states (10) belong at most to $`W`$. The range for the parameter $`p`$, in which $`𝒲_{W_2}`$ detects $`\rho `$, i.e. Tr$`(𝒲_{W_2}\rho )<0`$, is found to be $`3/5<p1`$, and is bigger than the one found by using $`𝒲_{W_1}`$. Taking any $`p`$ which has a finite distance to the border of this interval, i.e. $`p3/5>\mathrm{\Delta }`$ and $`1p>\mathrm{\Delta }`$, it is always possible to find a finite region around $`\rho `$ which still belongs to the $`WB`$–class. This can be seen by considering

$$\stackrel{~}{\rho }=(1ϵ)\left[\frac{1p}{8}1𝐥+pP_W\right]+ϵ\sigma ,$$
(11)
where $`\sigma `$ is an arbitrary density matrix, which covers all directions of possible deviations from $`\rho `$ in the operator space. In the worst case $`\sigma `$ is orthogonal to $`P_{GHZ}`$, so that Tr$`(P_{GHZ}\sigma )=0`$, and therefore Tr$`(𝒲_{W_2}\stackrel{~}{\rho })=(1ϵ)\text{Tr}(𝒲_{W_2}\rho )+ϵ/2`$. As long as the relation $`ϵ<(5p3)/(5p+1)`$ holds, the corresponding state $`\stackrel{~}{\rho }`$ is still detected by $`𝒲_{W_2}`$. Moreover, one can also find a finite $`ϵ^{}`$ such that if $`ϵ<ϵ^{}`$, then $`\stackrel{~}{\rho }`$ is in the $`W`$–class. The bound $`ϵ^{}`$ is obtained, for instance, by demanding that $`(1ϵ^{})(1p)1𝐥/8+ϵ^{}\sigma `$ is biseparable. The intersection of the two intervals gives a finite range for $`ϵ`$ where the state $`\stackrel{~}{\rho }`$ is in the $`WB`$–class. This proves that the set of mixed W $``$ B–states contains a ball, i.e. is not of measure zero.

We discuss now some possible consequences of our results for PPTES of three qubits, for which the partial transposes $`\rho ^{T_A}`$, $`\rho ^{T_B}`$ and $`\rho ^{T_C}`$ are positive. Any of these states can be decomposed as:

$$\rho =\lambda _S\rho _S+(1\lambda _S)\delta ,$$
(12)
where $`\rho _S`$ is a separable state and $`\delta `$ is an edge state. We conjecture that PPTES cannot belong to the GHZ $``$ W–class, i.e. they are at most in the $`W`$–class. This conjecture is rigorous for states that have edge states with low ranks in the above decomposition. It was shown in  that for bipartite systems in $`𝒞^2𝒞^N`$, the rank of PPTES must be larger than $`N`$, and if $`r(\rho )N`$ and $`\rho ^{T_A}0`$, then the state $`\rho `$ is separable. Thus, any PPTES of three–qubits with $`r(\rho )4`$ is biseparable with respect to any partition; an example of such states are the UPB–states from Ref. .

For the case of higher ranks we can only give some support for our conjecture. We proceed as in , and observe first that it suffices to prove the conjecture for the edge states. For these states, the sum of ranks satisfies $`r(\delta )+r(\delta ^{T_A})+r(\delta ^{T_B})+r(\delta ^{T_C})28`$. Any PPT entangled state can only be detected by a non-decomposable entanglement witness, which in the case of tripartite systems has the canonical form $`𝒲_{nd}=𝒲_dϵ1𝐥`$ where $`𝒲_d=P+Q_X^{T_X}`$ is a decomposable operator with $`P,Q_X0`$, $`R(P)=K(\delta )`$, $`R(Q_X)=K(\delta ^{T_X})`$ for some edge state $`\delta `$, and $`X=A,B,C`$ . We restrict ourselves to edge states with the maximal sum of ranks, i.e. states $`\delta `$ with $`(r(\delta ),r(\delta ^{T_A}),r(\delta ^{T_B}),r(\delta ^{T_C}))=(8,8,7,5),(8,8,6,6),(8,7,7,6),(7,7,7,7)`$ and permutations. Indeed, if the conjecture is true for these states, it will be true for all edge states, and thus for all PPTES, since the edge states with maximal sum of ranks are dense in the set of all edge states . We conjecture that for the case of edge states with maximal sum of ranks it is always possible to find a pure W–type vector, $`|\varphi _W`$, such that for any non-decomposable witness $`𝒲_{nd}`$ of $`\delta `$, $`\varphi _W|𝒲_d|\varphi _W0`$, so that $`\varphi _W|𝒲_{nd}|\varphi _W<0`$. That means $`𝒲_{nd}`$ cannot be a GHZ–witness, so the edge state $`\delta `$ belongs to the $`W`$–class. If this holds for any $`\delta `$ it implies that all PPTES belong to the $`W`$–class.
Any W–vector can be obtained by local invertible operations applied to $`|W`$ i.e. can be written as $`|\varphi _W=\alpha _A|e_2,f_1,g_1+\alpha _B|e_1,f_2,g_1+\alpha _C|e_1,f_1,g_2`$. We denote $`|\mathrm{\Phi }_A=|e_2^{},f_1,g_1`$, $`|\mathrm{\Psi }_A=\alpha _B|e_1^{},f_2,g_1+\alpha _C|e_1^{},f_1,g_2`$, $`|\mathrm{\Phi }_B=|e_1,f_2^{},g_1`$, $`|\mathrm{\Psi }_B=\alpha _A|e_2,f_1^{},g_1+\alpha _C|e_1,f_1^{},g_2`$, $`|\mathrm{\Phi }_C=|e_1,f_1,g_2^{}`$, $`|\mathrm{\Psi }_C=\alpha _A|e_2,f_1,g_1^{}+\alpha _B|e_1,f_2,g_1^{}`$. In order to fulfill the condition $`\varphi _W|𝒲_d|\varphi _W0`$ we demand that $`Q_X|\mathrm{\Phi }_X=0`$; $`P|\varphi _W=0`$, and $`Q_X|\mathrm{\Psi }_X=0`$ for $`X=A,B,C`$. The latter 4 conditions form 4 linear homogeneous equations for the $`\alpha _X`$’s, whose solutions exist if two $`3\times 3`$ determinants vanish. Together with the first 3 conditions this gives at most 5 equations in the case $`r(\delta )<8`$, and 6 equations in the worst case $`r(\delta )=8`$, for the 6 complex parameters characterizing $`|e_i,|f_i`$, and $`|g_i`$, with $`i=1,2`$. For $`r(\delta )<8`$ ($`r(\delta )=8`$) one expects here a one complex parameter (finite, but large) family of solutions. At the same time $`\varphi _W|𝒲_d|\varphi _W=2\mathrm{Re}_X\alpha _X\mathrm{\Psi }_X^_X|Q_X^{T_X}|\mathrm{\Phi }_X^_X`$, (where $`|\mathrm{\Phi }^_X`$ denotes partial complex conjugation with respect to $`X`$) i.e. is a hermitian form of $`\alpha _X`$’s, whose diagonal elements vanish, since $`|\mathrm{\Psi }_X`$ does not depend on $`\alpha _X`$. Employing the freedom of choosing the solutions from the family, one expects to find at least one with $`\varphi _W|𝒲_d|\varphi _W0`$. In this way we obtain the W–vector we were looking for. For the cases $`(6,8,8,6)`$ and $`(5,8,8,7)`$, a similar argument indeed shows that there should exist a biseparable state, $`|\psi _B`$, such that $`\psi _B|𝒲_{nd}|\psi _B<0`$. Note that the above method of searching $`|\psi _W`$ ($`|\psi _B`$) for a given $`\delta `$, if successful, provides a sufficient condition for $`\delta `$ to belong to the $`W`$–class ($`B`$–class).

Finally, we present an example for a PPTES entangled edge state with ranks (7,7,7,7). We introduce

$$\rho =\frac{1}{n}\left(\begin{array}{cccccccc}1& 0& 0& 0& 0& 0& 0& 1\\ 0& a& 0& 0& 0& 0& 0& 0\\ 0& 0& b& 0& 0& 0& 0& 0\\ 0& 0& 0& c& 0& 0& 0& 0\\ 0& 0& 0& 0& \frac{1}{c}& 0& 0& 0\\ 0& 0& 0& 0& 0& \frac{1}{b}& 0& 0\\ 0& 0& 0& 0& 0& 0& \frac{1}{a}& 0\\ 1& 0& 0& 0& 0& 0& 0& 1\end{array}\right)$$
(13)
with $`a,b,c>0`$ and $`n=2+a+1/a+b+1/b+c+1/c`$. The basis is $`\{000,001,010,011,100,101,110,111\}`$. This density matrix has a positive partial transpose with respect to each subsystem. One sees immediately that $`r(\rho )=r(\rho ^{T_A})=r(\rho ^{T_B})=r(\rho ^{T_{AB}})=7`$. In order to check that $`\rho `$ is a PPT entangled edge state, one has to prove that it is impossible to find a product vector $`|\varphi R(\rho )`$, such that at the same time $`|\varphi ^_XR(\rho ^{T_X})`$ for $`X=A,B,C`$. This, indeed, is not possible, as one readily concludes by looking at the kernels directly: one cannot find a product vector $`|\varphi `$ that is orthogonal to $`|\mathrm{\hspace{0.17em}000}|\mathrm{\hspace{0.17em}111}`$, whereas at the same time $`|\varphi ^_A|\mathrm{\hspace{0.17em}011}c|\mathrm{\hspace{0.17em}100}`$, $`|\varphi ^_B|\mathrm{\hspace{0.17em}010}b|\mathrm{\hspace{0.17em}101}`$, and $`|\varphi ^_C|\mathrm{\hspace{0.17em}001}a|\mathrm{\hspace{0.17em}110}`$, unless the condition $`ab=c`$ is fulfilled. Thus, for generic $`a,b,c`$ we have found a family of bound PPT entangled edge states of three qubits with maximal sum of ranks. By direct inspection we observe that $`\rho `$ fulfills our conjecture, and is biseparable with respect to any partition. It can be written e.g. as a sum of separable projectors and a B–state acting in the $`2\times 2`$ subspace spanned by Alice’s space and the vectors $`|\mathrm{\hspace{0.17em}00}`$ and $`|\mathrm{\hspace{0.17em}11}`$ in Bob’s–Charlie’s space.

To summarize, we show that the set of density matrices for three qubits has an “onion” structure (see Fig.1) and contains convex compact subsets of states belonging to the separable $`S`$, biseparable $`B`$, $`W`$– and $`GHZ`$–class, respectively. We provide the canonical way of constructing witness operators for the $`GHZ`$– and $`W`$–class, and give the first examples of such witnesses. The study of the family of tripartite states given in Eq. (10 ) allows us to prove that the $`W`$–class is not of measure zero. We conjecture and give some evidence that all PPTES of three–qubit systems do not require GHZ–type pure states for their formation. We formulate a sufficient condition which allows to check constructively if a state belongs to the $`W`$–class ($`B`$–class). Finally, we present a family of PPT entangled edge states of three qubits with maximal sum of ranks.

This work has been supported by DFG (SFB 407 and Schwerpunkt “Quanteninformationsverarbeitung”), the ESF-Programme PESC, and the EU IST-Programme EQUIP. AA thanks the University of Hannover for hospitality, E. Jané for useful comments and the Spanish MEC (AP-98) for financial support.


format: text

text:
# Choosing integration points for QCD calculations by numerical integration

## I Introduction

There is an important class of computer programs that do calculations in quantum chromodynamics (QCD) in which the calculation is performed at next-to-leading order in perturbation theory and allows for the determination of a variety of characteristics of the final state. This paper concerns a significant technical issue that arises in such programs when a “completely numerical” integration algorithm is used: how should one choose the integration points?

I consider the calculation of “three-jet-like” observables in $`e^+e^{}`$ annihilation. A given program in the class considered can be used to calculate a jet cross section (with any infrared safe choice of jet definition) or observables like the thrust distribution. Such a program generates random partonic events consisting of three or four final state quarks, antiquarks, and gluons. Each event comes with a calculated weight. A separate routine then calculates the contribution to the desired observable for each event, averaging over the events with their weights.

The weights are treated as probabilities. However, these weights can be both positive or negative. This is an almost inevitable consequence of quantum mechanics. The calculated observable is proportional to the square of a quantum amplitude and is thus positive. However, as soon as one divides the amplitude into pieces for purposes of calculation, one finds that while the square of each piece is positive, the interference terms between different pieces can have either sign. Thus the kind of program discussed here stands in contrast to the tree-level event generators in which, by simplifying the physics, one can generally arrange to have all the weights be positive, or, even, to be all equal to 1.

To understand the algorithms used in the class of programs described above, it is best to think of the programs as performing integrations over momenta in which the quantum matrix elements and the measurement functions form the integrand. There are two basic algorithms for performing the integrations. The older is due to Ellis, Ross, and Terrano (ERT) . (The first general purpose implementation of this method for three-jet-like observables in $`e^+e^{}`$ annihilation was that of Kunszt and Nason ). In this method, some of the integrations are performed analytically ahead of time. The other integrations are performed numerically by the Monte Carlo method. The integrations are divergent and are regulated by analytical continuation to $`32ϵ`$ space dimensions and a scheme of subtractions or cutoffs. The second method is much newer . In this method, all of the integrations are done by Monte Carlo numerical integration. With this method, the integrals are all convergent (after removal of the ultraviolet divergences by a straightforward renormalization procedure).

Since this method is quite new, one cannot yet say for what problems it might do better than the now standard ERT method. I can point out that the current incarnation of the numerical method has convergence problems for three jet quantities that receive important contributions from final states that are close to the two jet limit. For example, one gets good results for the three jet cross section or for the thrust distribution away from $`T=1`$, but poorly converging results for the energy-energy correlation function. On the other hand, the numerical method offers evident advantages in flexibility to modify the integrand. Thus one should be able to attack problems that are not accessible to the ERT method.

The numerical integration method exists as computer code  with accompanying technical notes  and many of the basic ideas behind it have been described in the two papers . The present paper is devoted to an exposition of the choice of integration points needed for this method. This subject may seem unimportant. Furthermore, it may seem to be a just a part of the numerical algorithms branch of computer science and thus to be uninteresting from the point of view of physics. However, practitioners of the art of numerical integrations in physics know that the choice of integration points is extremely important. If everything else is perfect but the integration points are badly chosen, the computer will, like Sisyphus, be faced with an unending task because the precision of the answer will not improve as the computational time increases. Furthermore, the problem is a problem in physics. All that the theory of numerical integration tells us is that the density of integration points should be well matched to the structure of the function to be integrated. The function to be integrated comes from Feynman diagrams, so we are really faced with the problem of understanding the structure of the quantum scattering processes. Perhaps surprisingly, the aspects of quantum scattering that are important in numerical method are completely different than the aspects that are important in the ERT method.

In Sec. II below, I briefly review the basics of the numerical integration method with the aim of setting the notation and making it possible to read this paper independently of the previous two papers. Then, in the sections that follow, I turn to the main issue of this paper, the choice of integration points. I try to keep this discussion succinct. The main features of the method are covered, but implementation details are left to the program code and its accompanying technical notes .

The exposition begins in Sec. III with a review of the general principles of Monte Carlo integration as they apply to this calculation. Sec. IV is the heart of the paper and covers the organization of the method used to choose the most integration points in the most important integration regions. Part of the method used considers the possible three parton final states. This part is explained in Sec. V. The other parts of the method concern the singularities of virtual graphs and their relation to the energy denominators in 2 parton to $`n`$ parton scattering. The various possibilities are covered in Secs. VI through IX. Of special importance is the elliptical coordinate system introduced in Sec. VI. A brief summary is presented in Sec. X.

## II Review of the numerical method

Let us begin with a precise statement of the problem. We consider an infrared safe three-jet-like observable in $`e^+e^{}\mathrm{ℎ𝑎𝑑𝑟𝑜𝑛𝑠}`$ such as a particular moment of the thrust distribution. The observable can be expanded in powers of $`\alpha _s/\pi `$,

$$\sigma =\underset{n}{}\sigma ^{[n]},\sigma ^{[n]}\left(\alpha _s/\pi \right)^n.$$
(1)
The order $`\alpha _s^2`$ contribution has the form

$`\sigma ^{[2]}`$ $`=`$ $`{\displaystyle \frac{1}{2!}}{\displaystyle 𝑑\stackrel{}{p}_1𝑑\stackrel{}{p}_2\frac{d\sigma _2^{[2]}}{d\stackrel{}{p}_1d\stackrel{}{p}_2}𝒮_2(\stackrel{}{p}_1,\stackrel{}{p}_2)}`$ (4)
$`+{\displaystyle \frac{1}{3!}}{\displaystyle 𝑑\stackrel{}{p}_1𝑑\stackrel{}{p}_2𝑑\stackrel{}{p}_3\frac{d\sigma _3^{[2]}}{d\stackrel{}{p}_1d\stackrel{}{p}_2d\stackrel{}{p}_3}𝒮_3(\stackrel{}{p}_1,\stackrel{}{p}_2,\stackrel{}{p}_3)}`$
$`+{\displaystyle \frac{1}{4!}}{\displaystyle 𝑑\stackrel{}{p}_1𝑑\stackrel{}{p}_2𝑑\stackrel{}{p}_3𝑑\stackrel{}{p}_4\frac{d\sigma _4^{[2]}}{d\stackrel{}{p}_1d\stackrel{}{p}_2d\stackrel{}{p}_3d\stackrel{}{p}_4}𝒮_4(\stackrel{}{p}_1,\stackrel{}{p}_2,\stackrel{}{p}_3,\stackrel{}{p}_4)}.`$
Here the $`d\sigma _n^{[2]}`$ are the order $`\alpha _s^2`$ contributions to the parton level cross section, calculated with zero quark masses. Each contains momentum and energy conserving delta functions. The $`d\sigma _n^{[2]}`$ include ultraviolet renormalization in the $`\overline{\mathrm{MS}}`$ scheme. The functions $`𝒮`$ describe the measurable quantity to be calculated. We wish to calculate a “three-jet-like” quantity. That is, $`𝒮_2=0`$. The normalization is such that $`𝒮_n=1`$ for $`n=2,3,4`$ would give the order $`\alpha _s^2`$ perturbative contribution the total cross section. There are, of course, infrared divergences associated with Eq. (4). For now, we may simply suppose that an infrared cutoff has been supplied.

The measurement, as specified by the functions $`𝒮_n`$, is to be infrared safe, as described in Ref. : the $`𝒮_n`$ are smooth functions of the parton momenta and

$$𝒮_{n+1}(\stackrel{}{p}_1,\mathrm{},\lambda \stackrel{}{p}_n,(1\lambda )\stackrel{}{p}_n)=𝒮_n(\stackrel{}{p}_1,\mathrm{},\stackrel{}{p}_n)$$
(5)
for $`0\lambda <1`$. That is, collinear splittings and soft particles do not affect the measurement.

It is convenient to calculate a quantity that is dimensionless. Let the functions $`𝒮_n`$ be dimensionless and eliminate the remaining dimensionality in the problem by dividing by $`\sigma _0`$, the total $`e^+e^{}`$ cross section at the Born level. Let us also remove the factor of $`(\alpha _s/\pi )^2`$. Thus, we calculate

$$=\frac{\sigma ^{[2]}}{\sigma _0(\alpha _s/\pi )^2}.$$
(6)
Let us now see how to set up the calculation of $``$ in a convenient form. We note that $``$ is a function of the c.m. energy $`\sqrt{s}`$ and the $`\overline{\mathrm{MS}}`$ renormalization scale $`\mu `$. We will choose $`\mu `$ to be proportional to $`\sqrt{s}`$: $`\mu =A_{UV}\sqrt{s}`$. Then $``$ depends on $`A`$. But, because it is dimensionless, it is independent of $`\sqrt{s}`$. This allows us to write

$$=_0^{\mathrm{}}𝑑\sqrt{s}h(\sqrt{s})(A_{UV},\sqrt{s}),$$
(7)
where $`h`$ is any function with

$$_0^{\mathrm{}}𝑑\sqrt{s}h(\sqrt{s})=1.$$
(8)
The quantity $``$ can be expressed in terms of cut Feynman diagrams, as in Fig. 1. The dots where the parton lines cross the cut represent the function $`𝒮_n(\stackrel{}{p}_1,\mathrm{},\stackrel{}{p}_n)`$. Each diagram is a three loop diagram, so we have integrations over loop momenta $`l_1^\mu `$, $`l_2^\mu `$ and $`l_3^\mu `$. We first perform the energy integrations. For the graphs in which four parton lines cross the cut, there are four mass-shell delta functions $`\delta (k_J^2)`$. These delta functions eliminate the three energy integrals over $`l_1^0`$, $`l_2^0`$, and $`l_3^0`$ as well as the integral (8) over $`\sqrt{s}`$. For the graphs in which three parton lines cross the cut, we can eliminate the integration over $`\sqrt{s}`$ and two of the $`l_J^0`$ integrals. One integral over the energy $`E`$ in the virtual loop remains. We perform this integration by closing the integration contour in the lower half $`E`$ plane. This gives a sum of terms obtained from the original integrand by some simple algebraic substitutions. Having performed the energy integrations, we are left with an integral of the form

$$=\underset{G}{}𝑑\stackrel{}{l}_1𝑑\stackrel{}{l}_2𝑑\stackrel{}{l}_3\underset{C}{}g(G,C;\stackrel{}{l}_1,\stackrel{}{l}_2,\stackrel{}{l}_3).$$
(9)
Here there is a sum over graphs $`G`$ (of which one is shown in Fig. 1) and there is a sum over the possible cuts $`C`$ of a given graph. The problem of calculating $``$ is now set up in a convenient form for calculation.

If we were using the Ellis-Ross-Terrano method, we would put the sum over cuts outside of the integrals in Eq. (9). For those cuts $`C`$ that have three partons in the final state, there is a virtual loop. We can arrange that one of the loop momenta, say $`\stackrel{}{l}_1`$, goes around this virtual loop. The essence of the ERT method is to perform the integration over the virtual loop momentum analytically ahead of time. The integration is often ultraviolet divergent, but the ultraviolet divergence is easily removed by a renormalization subtraction. The integration is also typically infrared divergent. This divergence is regulated by working in $`32ϵ`$ space dimensions and then taking $`ϵ0`$ while dropping the $`1/ϵ^n`$ contributions (after proving that they cancel against other contributions). After the $`\stackrel{}{l}_1`$ integration has been performed analytically, the integrations over $`\stackrel{}{l}_2`$ and $`\stackrel{}{l}_3`$ can be performed numerically. For the cuts $`C`$ that have four partons in the final state, there are also infrared divergences. One uses either a “phase space slicing” or a “subtraction” procedure to get rid of these divergences, cancelling the $`1/ϵ^n`$ pieces against the $`1/ϵ^n`$ pieces from the virtual graphs. In the end, we are left with an integral $`𝑑\stackrel{}{l}_1𝑑\stackrel{}{l}_2𝑑\stackrel{}{l}_3`$ in exactly three space dimensions that can be performed numerically.

In the numerical method, we keep the sum over cuts $`C`$ inside the integrations. We take care of the ultraviolet divergences by simple renormalization subtractions on the integrand. We make certain deformations on the integration contours so as to keep away from poles of the form $`1/[E_FE_I\pm iϵ]`$, where $`E_F`$ is the energy of the final state and $`E_I`$ is the energy of an intermediate state. Then the integrals are all convergent and we calculate them by Monte Carlo numerical integration.

Let us now look at the contour deformation in a little more detail. We denote the momenta $`\{\stackrel{}{l}_1,\stackrel{}{l}_2,\stackrel{}{l}_3\}`$ collectively by $`l`$ whenever we do not need a more detailed description. Thus

$$=\underset{G}{}𝑑l\underset{C}{}g(G,C;l).$$
(10)
For cuts $`C`$ that leave a virtual loop integration, there are singularities in the integrand of the form $`E_FE_I+iϵ`$ (or $`E_FE_Iiϵ`$ if the loop is in the complex conjugate amplitude to the right of the cut). Here $`E_F`$ is the energy of the final state defined by the cut $`C`$ and $`E_I`$ is the energy of a possible intermediate state. We will examine the nature of these “scattering singularities” in some detail in the next section. For now, all we need to know is that $`E_FE_I=0`$ on a surface in the space of $`\stackrel{}{l}_1`$ for fixed $`\stackrel{}{l}_2`$ and $`\stackrel{}{l}_3`$ if we pick $`\stackrel{}{l}_1`$ to be the momentum that flows around the virtual loop. These singularities do not create divergences. The Feynman rules provide us with the $`iϵ`$ prescriptions that tell us what to do about the singularities: we should deform the integration contour into the complex $`\stackrel{}{l}_1`$ space so as to keep away from them. Thus we write our integral in the form

$$=\underset{G}{}𝑑l\underset{C}{}𝒥(G,C;l)g(G,C;l+i\kappa (G,C;l)).$$
(11)
Here $`i\kappa `$ is a purely imaginary nine-dimensional vector that we add to the real nine-dimensional vector $`l`$ to make a complex nine-dimensional vector. The imaginary part $`\kappa `$ depends on the real part $`l`$, so that when we integrate over $`l`$, the complex vector $`l+i\kappa `$ lies on a surface, the integration contour, that is moved away from the real subspace. When we thus deform the contour, we supply a jacobian $`𝒥=det((l+i\kappa )/l)`$. (See Ref.  for details.)

The amount of deformation $`\kappa `$ depends on the graph $`G`$ and, more significantly, the cut $`C`$. For cuts $`C`$ that leave no virtual loop, each of the momenta $`\stackrel{}{l}_1`$, $`\stackrel{}{l}_2`$, and $`\stackrel{}{l}_3`$ flows through the final state. For practical reasons, we want the final state momenta to be real. Thus we set $`\kappa =0`$ for cuts $`C`$ that leave no virtual loop. On the other hand, when the cut $`C`$ does leave a virtual loop, we choose a non-zero $`\kappa `$. We must, however, be careful. When $`\kappa =0`$ there are singularities in $`g`$ on certain surfaces that correspond to collinear parton momenta. These singularities cancel between $`g`$ for one cut $`C`$ and $`g`$ for another. This cancellation would be destroyed if, for $`l`$ approaching the collinear singularity, $`\kappa =0`$ for one of these cuts but not for the other. For this reason, we insist that for all cuts $`C`$, $`\kappa 0`$ as $`l`$ approaches one of the collinear singularities. The details can be found in Ref. . All that is important here is that $`\kappa 0`$ quadratically with the distance to a collinear singularity.

Much has been left out in this brief overview, but we should now have enough background to see the requirements for a sensible choice of integration points.

## III General principles for the choice of points

In the following sections, we discuss the choice of integration points for the evaluation of a given graph. In this section, we summarize the general principles of Monte Carlo integration as they apply to this calculation.

We wish to perform an integral of the form

$$=\underset{G}{}𝑑lf(G;l)$$
(12)
where

$$f(G;l)=\mathrm{}\left\{\underset{C}{}𝒥(G,C;l)g(G,C;l+i\kappa (G,C;l))\right\}.$$
(13)
(We can take the real part because we know in advance that $``$ is real.) Using the Monte Carlo technique, for each graph $`G`$ we choose a large number $`\xi (G)N`$ of points $`l`$, with $`_G\xi (G)=1`$ so that the total number of points is $`N`$. We sample the points $`l`$ for graph $`G`$ at random with a density $`\rho (G;l)`$, with normalization $`𝑑l\rho (G;l)=1`$. Then the integral is approximated by

$$_N=\frac{1}{N}\underset{G}{}\frac{1}{\xi (G)}\underset{j=1}{\overset{\xi (G)N}{}}\frac{f(l_j)}{\rho (G;l_j)}.$$
(14)
This is an approximation for the integral in the sense that, if we repeat the procedure a lot of times, the expectation value for $`_N`$ is

$$_N=.$$
(15)
The expected r.m.s. error is $``$, where

$$^2=\left(_N\right)^2.$$
(16)
With a little manipulation, one can rewrite this as

$$^2=\frac{1}{N}\left\{\underset{G}{}\xi (G)\left[\frac{\mathrm{\Delta }(G)}{\xi (G)}\underset{L}{}\mathrm{\Delta }(L)\right]^2+\left(\underset{L}{}\mathrm{\Delta }(L)\right)^2\right\}.$$
(17)
where

$`\mathrm{\Delta }(G)^2`$ $`=`$ $`{\displaystyle 𝑑l\rho (G;l)\left[\frac{|f(G;l)|}{\rho (G;l)}𝑑l^{}|f(G;l^{})|\right]^2}`$ (19)
$`+{\displaystyle 𝑑l\left(|f(G;l)|+f(G;l)\right)\times 𝑑l^{}\left(|f(G;l^{})|f(G;l^{})\right)}.`$
We see, first of all, that the expected error decreases proportionally to $`1/\sqrt{N}`$. Second, we see that for a given choice of the density functions $`\rho (G;l)`$, the ideal choice of the factors $`\xi (G)`$ is $`\xi (G)\mathrm{\Delta }(G)`$. This is, in fact, easy to implement. Third, the ideal choice of $`\rho (G;l)`$ is $`\rho (G;l)|f(G;l)|`$. This is, in fact, impossible to implement.

Although it is not possible to choose $`\rho (G;l)|f(G;l)|`$, at least we can choose it so that $`|f(G;l)|/\rho (G;l)`$ is not singular at the singularities of $`|f(l)|`$. Furthermore, we can try to make $`\rho (G;l)`$ big at places where we know that $`|f(G;l)|`$ is big.

We will build the general sampling method out of elementary sampling methods. That is, we will find a number of simple methods to choose points $`l`$ for our graph. Let the density of points for the $`i`$th elementary sampling method be $`\rho _i(G;l)`$ (normalized to $`𝑑l\rho _i(G;l)=1`$). Then we devote a fraction $`\lambda _i(G)`$ of the points to the choice with density $`\rho _i(G;l)`$ and obtain a net density

$$\rho (G;l)=\underset{i}{}\lambda _i(G)\rho _i(G;l).$$
(20)
In this way, we make the sampling problem manageable. If we know that $`|f(G;l)|`$ is big near a certain point or surface in the space of loop momenta, we can design one of the elementary sampling methods so that the corresponding $`\rho _i(G;l)`$ is big there. In undertaking this task, we do not have to simultaneously arrange that $`\rho _i(G;l)`$ be big at other places where $`|f(G;l)|`$ is big.

In the following sections, we discuss the elementary sampling methods. We imagine that we are dealing with only one specific graph $`G`$, so we suppress the index $`G`$ in the notation.

## IV Organization of the sampling method

We consider a fixed (uncut) graph, dropping references to the label $`G`$ of the graph. As mentioned in the previous section, we sample points $`l`$ in the space of loop momenta according to several elementary sampling methods, each labelled by and index $`i`$ and having a density of points $`\rho _i(l)`$. The net density is then $`\rho (l)=\lambda _i\rho _i(l)`$.

We first address the identification of loop momenta. There are eight propagator momenta $`\stackrel{}{k}_P`$. (See, for example, Fig. 1.) The three loop momenta $`\stackrel{}{l}_L`$ can in general be any three linearly independent linear combinations of the $`\stackrel{}{k}_P`$. We will keep the choice simple by choosing three of the $`\stackrel{}{k}_P`$ to be the loop momenta. However, this still leaves us with the choice of which three of the $`\stackrel{}{k}_P`$ should be loop momenta. It proves convenient to make different choices for different elementary sampling methods. We specify the choice by specifying a triplet of integers $`\{Q(1),Q(2),Q(3)\}`$ such that $`\stackrel{}{l}_J`$ is $`\stackrel{}{k}_{Q(J)}`$. We call $`Q`$ an index set. Then the complete set of propagator matrices can related to the loop momenta by a matrix $`A`$:

$$\stackrel{}{k}_P=\underset{L=1}{\overset{3}{}}A_L^P\stackrel{}{l}_L.$$
(21)
Evidently, given the index set $`Q`$, the matrix $`A`$ can be constructed by using the topology of the graph.

Now we characterize certain surfaces, to be called scattering singularity surfaces, near which the integrand $`|f(l)|`$ is big. To do this, we consider the cuts of our graph in which three partons appear in the final state. For each such cut, there is a virtual loop. Let $`\stackrel{}{l}_1`$ be the loop momentum. More precisely, $`\stackrel{}{l}_1`$ will be the momentum, $`\stackrel{}{k}_{Q(1)}`$, of one of the propagators in the loop, but we defer for a moment specifying which one. As specified in Sec. II, the integration contour for $`\stackrel{}{l}_1`$ is deformed into the complex plane,<sup>*</sup><sup>*</sup>*We keep the momenta that enter the final state real. so that $`\stackrel{}{l}_1\stackrel{}{l}_{1,c}=\stackrel{}{l}_1+i\stackrel{}{\kappa }`$.

Before deformation, the integrand is singular on certain surfaces associated with simple scattering processes, the scattering singularity surfaces. How can this happen? There are four cases. Each case is illustrated by a Feynman graph in Fig. 2. The corresponding scattering singularity surface is illustrated in Fig. 3. The cases are

1. 2 to 2 (s). A virtual parton with momentum $`\stackrel{}{l}_1`$ merges with a virtual parton with momentum $`\stackrel{}{l}_2+\stackrel{}{l}_3\stackrel{}{l}_1`$ to produce a virtual parton with momentum $`\stackrel{}{l}_2+\stackrel{}{l}_3`$. This virtual parton divides into two partons with momenta $`\stackrel{}{l}_2`$ and $`\stackrel{}{l}_3`$ that enter the final state. In old-fashioned perturbation theory, there is an energy denominator $`|\stackrel{}{l}_2|+|\stackrel{}{l}_3||\stackrel{}{l}_1||\stackrel{}{l}_2+\stackrel{}{l}_3\stackrel{}{l}_1|+iϵ`$, which vanishes on the ellipsoid $`|\stackrel{}{l}_1|+|\stackrel{}{l}_2+\stackrel{}{l}_3\stackrel{}{l}_1|=|\stackrel{}{l}_2|+|\stackrel{}{l}_3|`$ in the space of $`\stackrel{}{l}_1`$. This is the scattering singularity surface. The contour deformation is non-zero on the entire scattering singularity surface.
2. 2 to 2 (t). A virtual parton with momentum $`\stackrel{}{l}_2+\stackrel{}{l}_1`$ scatters from a virtual parton with momentum $`\stackrel{}{l}_3\stackrel{}{l}_1`$ by exchanging a parton with momentum $`\stackrel{}{l}_1`$. Partons with momentum $`\stackrel{}{l}_2`$ and $`\stackrel{}{l}_3`$ emerge into the final state. There is a scattering singularity surface $`|\stackrel{}{l}_2+\stackrel{}{l}_1|+|\stackrel{}{l}_3\stackrel{}{l}_1|=|\stackrel{}{l}_2|+|\stackrel{}{l}_3|`$. In this case, there is also a singularity at $`\stackrel{}{l}_1=0`$ that arises from the propagator of the exchanged parton. This soft exchange singularity lies on the scattering singularity surface. Furthermore, in our treatment, the contour deformation vanishes at the soft exchange singularity so that, even after contour deformation, the integrand is singular there. This is, however, an integrable singularity.
3. 2 to 3. A virtual parton with momentum $`\stackrel{}{l}_1`$ collides with a virtual parton with momentum $`\stackrel{}{l}_1`$ to produce a final state with partons carrying momenta $`\stackrel{}{l}_2`$, $`\stackrel{}{l}_3`$, and $`\stackrel{}{l}_2\stackrel{}{l}_3`$. There is an energy denominator $`|\stackrel{}{l}_2|+|\stackrel{}{l}_3|+|\stackrel{}{l}_2+\stackrel{}{l}_3|2|\stackrel{}{l}_1|+iϵ`$, which is singular on the sphere $`|\stackrel{}{l}_1|=[|\stackrel{}{l}_2|+|\stackrel{}{l}_3|+|\stackrel{}{l}_2+\stackrel{}{l}_3|]/2`$. The contour deformation is non-zero on the entire scattering singularity surface.
4. 2 to 1. A virtual parton with momentum $`\stackrel{}{l}_1`$ combines with a virtual parton with momentum $`\stackrel{}{l}_2\stackrel{}{l}_1`$ to produce an on shell parton with momentum $`\stackrel{}{l}_2`$ that enters the final state. There is an energy denominator $`|\stackrel{}{l}_2||\stackrel{}{l}_1||\stackrel{}{l}_2\stackrel{}{l}_1|+iϵ`$, which vanishes on the line $`\stackrel{}{l}_1=x\stackrel{}{l}_2`$ with $`0x1`$. The contour deformation is chosen to vanish at this collinear singularity. As discussed in Ref. , the collinear singularity cancels when one sums over cuts. However, singularities at $`\stackrel{}{l}_1=0`$ an $`\stackrel{}{l}_1=\stackrel{}{l}_2`$ remain.

Thus for each cut there are a number of scattering singularity surfaces. There is a contour deformation that keeps the integrand from being singular except at special points on these surfaces. However, the integrand is sometimes still quite large near these surfaces. For this reason, we will design an elementary sampling method for each surface such that the density of points is big on the whole surface and singular at the special point if necessary.

There are two kinds of singularities associated with points $`\stackrel{}{k}_P=0`$ where a propagator momentum vanishes. First, there is a $`1/|\stackrel{}{k}_P|^2`$ singularity when the momentum of the exchanged parton in a 2 to 2 (t) scattering vanishes. The elementary sampling method associated with the 2 to 2 (t) scattering will be designed to take care of this kind of singularity. Second, there is a $`1/|\stackrel{}{k}_P|`$ singularity when any propagator momentum $`\stackrel{}{k}_P`$ approaches zero. This weak singularity arises simply because massless Lorentz invariant phase space is $`d\stackrel{}{k}/|\stackrel{}{k}|`$. As it turns out, these singularities in the density will automatically be provided by the combined sampling methods without having to specifically arrange for them.

The 2 to 1 scattering singularity surface is exceptional in the list above in that there is no singularity except for the two singular points at $`\stackrel{}{l}_1=0`$ and $`\stackrel{}{l}_2\stackrel{}{l}_1=0`$. These two singular points are of one or the other of two types mentioned above. Typically a 2 to 1 scattering subdiagram is part of a 2 to 2 (s) or 2 to 2 (t) scattering subdiagram and the two singularities are provided for by the 2 to 2 (s) or 2 to 2 (t) sampling methods. The exception is in the case of a self-energy subgraph that is connected to a final state parton. In this case, the 2 to 2 (t), 2 to 2 (s), and 2 to 3 sampling methods do not apply and we need an explicit 2 to 1 sampling method. Thus we will apply a 2 to 1 sampling method only in the case of a self-energy subgraph connected to a final state parton.

The previous argument indicates that for each scattering singularity surface in the space of the loop momentum $`\stackrel{}{l}_1`$ in a virtual loop, we should associate a method for choosing $`\stackrel{}{l}_1`$ that puts a high density of points near this surface. It is then useful to choose the other two loop momenta to be the momenta of two of the three partons in the final state. Thus the momenta of the final state partons are $`\stackrel{}{l}_2`$, $`\stackrel{}{l}_3`$, and $`\stackrel{}{l}_2\stackrel{}{l}_3`$. The integrand will be singular whenever the three partons approach a two jet configuration. Thus we choose the points $`\{\stackrel{}{l}_2,\stackrel{}{l}_3\}`$ so that the density of points is appropriately singular at the two jet configurations.

Thus we have a general scheme for organizing the elementary sampling methods. First, find all of the possible three parton cuts of the graph in question. Then, for each such final state cut, enumerate the scattering singularity surfaces that can occur, counting the 2 to 1 case only when the virtual loop is a self-energy diagram connected to a final state parton. There are five basic situations, as illustrated in Fig. 4. Each combination of a final state cut and a scattering singularity surface will correspond to an elementary sampling method.

One point should be emphasized for clarity. In the numerical method discussed in this paper, for each integration point, we compute $`f/\rho `$ where $`f`$ is the integrand and $`\rho `$ is the density of integration points. Contributions from all cuts of a given graph are calculated and summed to form the integrand $`f`$. The density $`\rho `$ is also a sum, with terms corresponding to each of the possible cuts of the graph. Thus there are independent sums over possible cuts in $`f`$ and in $`\rho `$.

Consider the five basic amplitudes with a virtual loop that are illustrated in Fig. 4. In two of these cases, there are more than one possible intermediate state involving the virtual loop and thus more than one scattering singularity surface. In these cases, it is important to understand how the scattering singularity surfaces fit together.

Let us examine, then, case $`(d)`$ in Fig. 4. Label the momenta as in Fig. 5. In Fig. 6 we see the scattering singularity surfaces in the space of the loop momentum $`\stackrel{}{l}\stackrel{}{l}_1`$ for an arbitrary fixed choice of the final state momenta $`\stackrel{}{p}_j`$. (We draw the figure for $`\stackrel{}{l}`$ in the plane of the $`\stackrel{}{p}_j`$). The vectors $`\stackrel{}{p}_1`$, $`\stackrel{}{p}_2`$, and $`\stackrel{}{p}_3`$ are indicated. The integration contour is deformed everywhere except along the line $`\stackrel{}{l}=x\stackrel{}{p}_3`$ with $`0<x<1`$, which is indicated as a solid line. There is an ellipsoidal 2 to 2 (s) surface, $`|\stackrel{}{l}|+|\stackrel{}{p}_1+\stackrel{}{p}_2\stackrel{}{l}|=|\stackrel{}{p}_1|+|\stackrel{}{p}_2|`$. There is also a spherical 3 to 2 surface, $`2|\stackrel{}{l}|=|\stackrel{}{p}_1|+|\stackrel{}{p}_2|+|\stackrel{}{p}_3|`$. The integration contour is deformed everywhere on both of these surfaces, so that the integrand is not singular anywhere. On the other hand, the deformation vanishes on the solid line, which can be rather near the ellipsoidal surface in the case that the angle between $`\stackrel{}{p}_1`$ and $`\stackrel{}{p}_2`$ is small. Thus the integrand can be rather big on the ellipsoid in this circumstance. Furthermore, the size of the integrand is enhanced where the ellipsoid is tangent to the spherical singularity surface. We will need to put an extra density of points in the region of this point of tangency.

Case $`(e)`$ in Fig. 4 just a little more complicated. Label the momenta as in Fig. 7. Fig. 6 shows the scattering singularity surfaces in the space of the loop momentum $`\stackrel{}{l}`$ for fixed final state momenta $`\stackrel{}{p}_j`$. The vectors $`\stackrel{}{p}_1`$, $`\stackrel{}{p}_2`$, and $`\stackrel{}{p}_3`$ are indicated. The integration contour is deformed everywhere except along the lines $`\stackrel{}{l}=x\stackrel{}{p}_1`$, $`\stackrel{}{l}\stackrel{}{p}_1=x\stackrel{}{p}_2`$ and $`\stackrel{}{l}\stackrel{}{p}_1\stackrel{}{p}_2=x\stackrel{}{p}_3`$, with $`0<x<1`$ in each case. These lines are are indicated as solid lines that form a triangle. There are two ellipsoidal 2 to 2 (t) surfaces, $`|\stackrel{}{l}|+|\stackrel{}{p}_1+\stackrel{}{p}_2\stackrel{}{l}|=|\stackrel{}{p}_1|+|\stackrel{}{p}_2|`$ and $`|\stackrel{}{l}\stackrel{}{p}_1|+|\stackrel{}{l}|=|\stackrel{}{p}_2|+|\stackrel{}{p}_3|`$. There is also a spherical 3 to 2 surface, $`2|\stackrel{}{l}|=|\stackrel{}{p}_1|+|\stackrel{}{p}_2|+|\stackrel{}{p}_3|`$. As in the previous case, $`(d)`$, we need an especially high density of integration points near where an ellipsoid is tangent to the sphere in the case that this point is near to the triangle, where the deformation vanishes. We have also the new feature that each of the ellipsoids shares a point with the triangle. At this point the contour deformation is zero, so the integrand is actually singular. This is the point where the momentum of the exchanged parton in the associated 2 to 2 (t) scattering vanishes. The density of integration points will have to have a corresponding singularity at these points.

This completes the description of the general organization of the sampling scheme. In the following sections, we outline its component parts.

## V Sampling for the final state

As discussed in the previous section, for each final state cut with three partons in the final state, we will define a number of elementary sampling methods. Part of each method is to choose with an appropriate density the momentum $`\stackrel{}{l}_1`$ that circles through the virtual loop that occurs when there are three final state partons. The other part is to choose the momenta $`\stackrel{}{l}_2`$, $`\stackrel{}{l}_3`$, and $`\stackrel{}{l}_2\stackrel{}{l}_3`$ of the final state partons. We address the choice of the final state momenta in this section.

Let the final state parton momenta be $`\stackrel{}{p}_1`$, $`\stackrel{}{p}_2`$, and $`\stackrel{}{p}_3`$, with $`_i\stackrel{}{p}_i=0`$. (Two of these will then be $`\stackrel{}{l}_2`$ and $`\stackrel{}{l}_3`$, but it doesn’t matter which ones.) Define $`E_{\mathrm{max}}={\scriptscriptstyle \frac{1}{2}}_i|\stackrel{}{p}_i|`$ and $`x_i=|\stackrel{}{p}_i|/E_{\mathrm{max}}`$. Then $`_ix_i=2`$ and $`0x_i1`$. We can write the integration over final state momenta as follows, using $`p_3^2=p_1^2+p_2^2+2p_1p_2\mathrm{cos}\theta _{12}`$ so that $`d\mathrm{cos}\theta _{12}=p_3dp_3/(p_1p_2)`$:

$``$ $`=`$ $`{\displaystyle 𝑑\stackrel{}{p}_1𝑑\stackrel{}{p}_2𝑑\stackrel{}{p}_3\delta (\stackrel{}{p}_i)f(p)}`$ (22)
$`=`$ $`{\displaystyle 2E_{\mathrm{max}}^6d\mathrm{ln}E_{\mathrm{max}}d\mathrm{cos}\theta _1d\varphi _1d\varphi _{12}x_1x_2x_3dx_1dx_2f(p)}.`$ (23)
Here $`\{\theta _1,\varphi _1\}`$ are the angles of $`\stackrel{}{p}_1`$ and $`\varphi _{12}`$ is the angle between $`\stackrel{}{p}_1\times \stackrel{}{p}_2`$ and $`\stackrel{}{p}_1\times \stackrel{}{e}`$, where $`\stackrel{}{e}`$ is an arbitrary reference vector.

If we sample $`\{\mathrm{ln}E_{\mathrm{max}},\mathrm{cos}\theta _1,\varphi _1,\varphi _{12}\}`$ with a density

$$\rho _A\frac{dN}{d\mathrm{ln}E_{\mathrm{max}}d\mathrm{cos}\theta _1d\varphi _1d\varphi _{12}}$$
(24)
and we sample $`\{x_1,x_2\}`$ with a density

$$\rho _B\frac{dN}{dx_1dx_2},$$
(25)
then the net density is

$$\rho \frac{dN}{d\stackrel{}{p}_1d\stackrel{}{p}_2}=\frac{\rho _A\rho _B}{2E_{\mathrm{max}}^6x_1x_2x_3}.$$
(26)
On symmetry grounds, $`\rho _A`$ should be independent of the angles. Its dependence on $`E_{\mathrm{max}}`$ is not too important. A convenient choice, normalized so that $`𝑑N=1`$, is

$$\frac{1}{\rho _A}=\frac{8\pi ^2}{3}\left(E_0/E_{\mathrm{max}}\right)^3\left[1+\left(E_{\mathrm{max}}/E_0\right)^3\right]^2,$$
(27)
where $`E_0`$ is a fixed parameter with dimensions of mass.

What should be the properties of $`\rho _B`$? We can expect the integrand $`f(p)`$ to be singular in any two jet region: when any of the $`\stackrel{}{p}_j`$ vanishes ($`x_j0`$) or when two of them become collinear (if $`\stackrel{}{p}_1`$ and $`\stackrel{}{p}_2`$ are nearly collinear, then $`x_31`$). The factor $`1/x_1x_2x_3`$ in Eq. (26) provides a weak singularity in the $`x_i0`$ regions. We can build in a weak singularity in the collinear regions by choosing

$$\rho _B=\frac{𝒩_x\underset{j}{\mathrm{max}}\sqrt{1x_j}}{_j\sqrt{1x_j}}.$$
(28)
The $`x_i0`$ regions are at the intersections of $`x_l1`$ regions, so these collinear singularities also enhance the soft singularity of the density. The normalization constant is fixed so that $`𝑑N=1`$ and is $`𝒩_x=1/(6\sqrt{2}\mathrm{arcsin}(1/\sqrt{3}))`$. One can, of course, change $`1/\sqrt{1x}`$ to $`1/(1x)^A`$ for any A smaller than 1.

It is a simple matter to actually choose points with the density $`\rho `$ in Eq. (26). Readers interested in the implementation details may consult the notes  that accompany the program code.

I can comment here about the singularities of the integrand, $`f`$. If, for the moment, we set the measurement function to 1, then the nominal behavior of $`f`$ is $`f1/x_j^3`$ for $`x_j0`$ and $`f1/(1x_j)`$ for $`x_j1`$. This nominal behavior is what one gets after performing the virtual loop integral, $`𝑑l_1`$, inside the $`\{\stackrel{}{p}_1,\stackrel{}{p}_2\}`$ integral . These singularities make the integration logarithmically divergent. Then the measurement function that is included in $`f`$ is required to vanish in the two jet limit, so that the integration becomes convergent. Now, in fact, if we look at the singularities in the $`\{\stackrel{}{p}_1,\stackrel{}{p}_2\}`$ integral before performing the virtual loop integral, the singularities are worse than this nominal behavior. For this reason, the suppression of the two jet region coming from the measurement function must be suitably strong in order to obtain good convergence. It remains for future research to arrange the calculation in such a way that the nominal behavior of the integrand $`f`$ as a function of $`\{\stackrel{}{p}_1,\stackrel{}{p}_2\}`$ is obtained.

## VI Sampling for 2 to 2 (s) scattering

In this section, we consider the sampling method associated with 2 to 2 (s) scattering. We need to choose points $`\stackrel{}{l}\stackrel{}{l}_1`$ appropriate to the following case:

> A virtual parton with momentum $`\stackrel{}{l}`$ merges with a virtual parton with momentum $`\stackrel{}{l}_2+\stackrel{}{l}_3\stackrel{}{l}`$ to produce a virtual parton with momentum $`\stackrel{}{l}_2+\stackrel{}{l}_3`$. This virtual parton divides into two partons with momenta $`\stackrel{}{l}_2`$ and $`\stackrel{}{l}_3`$ that enter the final state.

In this case, the scattering surface is the ellipsoid $`|\stackrel{}{l}|+|\stackrel{}{l}_2+\stackrel{}{l}_3\stackrel{}{l}|=|\stackrel{}{l}_2|+|\stackrel{}{l}_3|`$. The density of points should be large on this surface. In order to accomplish this, we use elliptical coordinates.

### A Elliptical coordinates

The elliptical coordinates $`A_+,A_{},\varphi `$ are defined as follows. First, define $`\kappa `$ by

$$2\kappa =|\stackrel{}{l}_2+\stackrel{}{l}_3|.$$
(29)
Now, define coordinates $`A_\pm `$ by

$$A_\pm =\frac{1}{2\kappa }\left(|\stackrel{}{l}|\pm |\stackrel{}{l}\stackrel{}{l}_2\stackrel{}{l}_3|\right).$$
(30)
Then

$$1<A_+,1<A_{}<1.$$
(31)
The coordinate $`A_+`$ is constant on elliptical surfaces, while $`A_{}`$ is constant on orthogonal hyperbolic surfaces. The scattering singularity is the ellipsoid

$$A_+=S_+,$$
(32)
where

$$S_\pm =\frac{1}{2\kappa }\left(|\stackrel{}{l}_2|\pm |\stackrel{}{l}_3|\right).$$
(33)
We need one more coordinate. Let $`\varphi `$ be the azimuthal angle of $`\stackrel{}{l}`$ in a coordinate system in which the $`z`$ axis lies in the direction of $`\stackrel{}{l}_2+\stackrel{}{l}_3`$ and $`\stackrel{}{l}_2`$ has zero $`y`$ component and a positive $`x`$ component. To state this precisely, define unit vectors

$`\stackrel{}{n}_z`$ $`=`$ $`{\displaystyle \frac{\stackrel{}{l}_2+\stackrel{}{l}_3}{|\stackrel{}{l}_2+\stackrel{}{l}_3|}}`$ (34)
$`\stackrel{}{n}_y`$ $`=`$ $`{\displaystyle \frac{\stackrel{}{l}_3\times \stackrel{}{l}_2}{|\stackrel{}{l}_3\times \stackrel{}{l}_2|}}`$ (35)
$`\stackrel{}{n}_x`$ $`=`$ $`\stackrel{}{n}_y\times \stackrel{}{n}_z.`$ (36)
Then

$$\varphi =\mathrm{arctan}\left(\stackrel{}{l}\stackrel{}{n}_y/\stackrel{}{l}\stackrel{}{n}_x\right).$$
(37)
The transformation from $`A_+,A_{},\varphi `$ to $`\stackrel{}{l}`$ is

$$\stackrel{}{l}=\frac{1}{2}(\stackrel{}{l}_2+\stackrel{}{l}_3)+l_T\mathrm{cos}\varphi \stackrel{}{n}_x+l_T\mathrm{sin}\varphi \stackrel{}{n}_y+z\stackrel{}{n}_z,$$
(38)
where

$`l_T`$ $`=`$ $`\kappa \left(A_+^21\right)^{1/2}\left(1A_{}^2\right)^{1/2}`$ (39)
$`z`$ $`=`$ $`\kappa A_+A_{}.`$ (40)
A straightforward calculation shows that the jacobian of the transformation is given by

$$d\stackrel{}{l}=\frac{dA_+dA_{}d\varphi }{\rho _{AA\varphi }},$$
(41)
where

$$\frac{1}{\rho _{AA\varphi }}=\kappa ^3(A_+^2A_{}^2).$$
(42)
The factor $`(A_+^2A_{}^2)`$ is convenient. It provides weak singularities in the density of points when $`A_+\pm A_{}0`$, which corresponds to $`\stackrel{}{l}0`$ and $`\stackrel{}{l}\stackrel{}{l}_2\stackrel{}{l}_30`$.

If we sample points in the variables $`\{A_+,A_{},\varphi \}`$ with a density

$$\rho ^{}=\frac{dN}{dA_+dA_{}d\varphi },$$
(43)
then the total density of points will be

$$\rho =\rho ^{}\times \rho _{AA\varphi },$$
(44)
where $`\rho _{AA\varphi }`$ is given in Eq. (42). We next address the question of what to choose for $`\rho ^{}`$.

### B The choice of $`A_+`$, $`A_{}`$ and $`\varphi `$

With what density

$$\rho ^{}=\frac{dN}{dA_+dA_{}d\varphi }$$
(45)
should we sample points in the variables $`\{A_+,A_{},\varphi \}`$? There is no unique answer, but consider a density function of the form

$$\rho ^{}=𝒩\frac{1}{1A_{}^2+\tau }\frac{1}{A_+(|A_+S_+|+\lambda )}.$$
(46)
Here $`\tau `$ and $`\lambda `$ are parameters to be specified and the normalization $`𝒩`$ is given by

$$\frac{1}{𝒩}=\frac{2\pi }{\sqrt{1+\tau }}\mathrm{ln}\left(\frac{\sqrt{1+\tau }+1}{\sqrt{1+\tau }1}\right)\left\{\frac{1}{S_++\lambda }\mathrm{ln}\left(\frac{S_+[S_+1+\lambda ]}{\lambda }\right)+\frac{1}{S_+\lambda }\mathrm{ln}\left(\frac{S_+}{\lambda }\right)\right\}.$$
(47)
With this ansatz, there is a high density of points near the scattering singularity surface, $`A_+=S_+`$. The width of the peak is $`\lambda `$. Thus $`\lambda `$ should be matched to the width of the peak in the integrand. If the singularity surface is far from the line $`A_+=1`$, the contour deformation is substantial and peak is broad. On the other hand, if the singularity surface is near to the line $`A_+=1`$, the contour deformation is small and peak is narrow. I estimate that the width of the peak in the integrand is of order $`(S_+1)^2`$ for $`S_+`$ near 1. For large $`S_+`$ it seems reasonable to choose a width of order $`S_+`$. Thus I take

$$\lambda =\frac{(S_+1)^2}{S_+}.$$
(48)
With the function $`\rho ^{}`$ in Eq. (46), there is an enhancement of the density of points near $`A_{}=\pm 1`$. This enhancement is built into the density in order to provide an extra density of points near where the ellipsoid in Fig. 6 is tangent to the sphere. As explained in Sec. IV, we need a high density of points here when $`\stackrel{}{l}_2`$ and $`\stackrel{}{l}_3`$ are nearly collinear, that is when $`S_+`$ is close to 1. I arrange for this by taking the parameter $`\tau `$ in Eq. (46) to be

$$\tau =\frac{S_+1}{S_+}.$$
(49)
It is easy to sample points $`\{A_+,A_{},\varphi \}`$ with the density $`\rho ^{}`$ in Eq. (46) as explained in .

## VII Sampling for 2 to 2 (t) scattering

In this section, we consider the sampling method associated with 2 to 2 (t) scattering. We need to choose points $`\stackrel{}{l}\stackrel{}{l}_1`$ appropriate to the following case:

> A virtual parton with momentum $`\stackrel{}{l}_2+\stackrel{}{l}`$ scatters from a virtual parton with momentum $`\stackrel{}{l}_3\stackrel{}{l}`$ by exchanging a parton with momentum $`\stackrel{}{l}`$. Partons with momentum $`\stackrel{}{l}_2`$ and $`\stackrel{}{l}_3`$ emerge into the final state.

In this case, there is a scattering singularity surface $`|\stackrel{}{l}_2+\stackrel{}{l}|+|\stackrel{}{l}_3\stackrel{}{l}|=|\stackrel{}{l}_2|+|\stackrel{}{l}_3|`$. There is a singularity on this surface at the point where the momentum $`\stackrel{}{l}`$ of the exchanged parton vanishes. The density of points should have a matching singularity at this point.

### A Elliptical coordinates

As in the previous section, we use elliptical coordinates $`A_+,A_{},\varphi `$. The role of $`\stackrel{}{l}`$ in the previous section is now played by $`\stackrel{}{l}+\stackrel{}{l}_2`$, so the formulas are a little different. We again define $`\kappa `$ by $`2\kappa =|\stackrel{}{l}_2+\stackrel{}{l}_3|.`$ Now, define coordinates $`A_\pm `$ by

$$A_\pm =\frac{1}{2\kappa }\left(|\stackrel{}{l}_2+\stackrel{}{l}|\pm |\stackrel{}{l}_3\stackrel{}{l}|\right).$$
(50)
Then $`1<A_+`$ and $`1<A_{}<1.`$ Define

$$S_\pm =\frac{1}{2\kappa }\left(|\stackrel{}{l}_2|\pm |\stackrel{}{l}_3|\right).$$
(51)
The scattering singularity is the ellipsoid

$$A_+=S_+$$
(52)
while the soft exchange singularity is at

$$\{A_+,A_{}\}=\{S_+,S_{}\}.$$
(53)
We need one more coordinate, an azimuthal angle $`\varphi `$. We define unit vectors $`\{\stackrel{}{n}_x,\stackrel{}{n}_y,\stackrel{}{n}_z\}`$ according to Eq. (36). Then we define

$$\varphi =\mathrm{arctan}\left((l_2+\stackrel{}{l})\stackrel{}{n}_y/(l_2+\stackrel{}{l})\stackrel{}{n}_x\right).$$
(54)
The jacobian of the transformation from $`\stackrel{}{l}`$ to $`\{A_+,A_{},\varphi \}`$ is again $`1/\rho _{AA\varphi }`$ where $`\rho _{AA\varphi }`$ is given in Eq. (42). If we sample points in the parameters $`\{A_+,A_{},\varphi \}`$ with a density $`\rho ^{}`$ then the total density of points will be $`\rho =\rho ^{}\times \rho _{AA\varphi }.`$ We next address the question of what to choose for $`\rho ^{}`$.

### B The choice of $`A_+`$, $`A_{}`$ and $`\varphi `$

In this subsection we specify a choice for the density

$$\rho ^{}=\frac{dN}{dA_+dA_{}d\varphi }$$
(55)
with which we sample in the variables $`\{A_+,A_{},\varphi \}`$. We begin by recognizing that we face two challenges. First, we would like to have a concentration of points near the surface $`A_+=S_+`$ with an especially high density near $`A_{}=\pm 1`$ if $`S_+`$ is near 1, as discussed with respect to the sampling for 2 to 2 (s). The second challenge is to make $`\rho ^{}`$ appropriately singular at $`\{A_+,A_{},\varphi \}=\{S_+,S_{},0\}`$. We thus take $`\rho ^{}`$ to have two parts

$$\rho ^{}=\alpha _{2t}\rho _s+(1\alpha _{2t})\rho _t,$$
(56)
where $`\alpha _{2t}`$ is a fixed parameter with $`0<\alpha _{2t}<1`$. We take $`\rho _s`$ to be given by the density (46) for the sampling for 2 to 2 (s). Then we have met the first challenge. It remains to design $`\rho _t`$ to address the second challenge.

The main idea is that there is a scattering singularity surface at $`A_+=S_+`$, so that the integrand has a factor

$$\frac{1}{A_+S_++i\eta }$$
(57)
where $`\eta `$ is the amount of deformation of the $`S_+`$ contour. The amount of contour deformation vanishes at the soft exchange singularity at $`\{A_+,A_{},\varphi \}=\{S_+,S_{},0\}`$ and, near to this point, $`\eta `$ is proportional to the square of the distance to $`\{S_+,S_{},0\}`$. Thus, on the surface $`A_+=S_+`$ near to the soft exchange singularity, we can estimate $`\eta `$ as $`\omega ^2`$ where

$$\omega |A_{}S_{}|+|\varphi |/\pi .$$
(58)
There is another factor of $`1/\omega `$ that arises from the propagator of the exchanged parton after we take real-virtual cancellations into account, as explained in detail in Ref. . Thus the absolute value of the integrand has a factor that can be estimated as

$$\frac{1}{\omega |A_+S_+|}$$
(59)
for $`\omega 1`$ and $`\omega ^2|A_+S_+|\omega `$. We want $`\rho _t`$ to have a singularity of the same nature, so that the integrand divided by the density of points is singularity free. Thus we take

$$\rho _t=\frac{𝒜}{[|A_+S_+|+\beta _{2t}S_+\omega ^2][|A_+S_+|+\beta _{2t}\gamma _{2t}S_+\omega ]}.$$
(60)
Here $`\beta _{2t}`$ and $`\gamma _{2t}`$ are fixed parameters and $`𝒜`$ is a rather complicated function of $`A_{}`$ and $`\varphi `$ (described below) that is of secondary importance. The main point is that $`\rho _t`$ behaves like (59), with cutoffs when $`|A_+S_+|`$ gets to be smaller than $`\omega ^2`$ or larger than $`\omega `$.

We have thus given $`\rho _t`$ an actual singularity at the point $`\{A_+,A_{},\varphi \}=\{S_+,S_{},0\}`$. This singularity in $`\rho _t`$ has a special structure such that matches the structure of the integrand $`f`$, so that $`f/\rho _t`$ is finite in the neighborhood of the soft exchange singularity. This point, which is important for the convergence of the numerical integration, was described in some detail in Ref. . In that paper, however, the construction was not as nice as one would like because the ellipsoidal coordinates were not used and the “ridgeline” of $`\rho `$ was placed on the tangent plane to the ellipsoid rather than on the ellipsoid itself.

To sample points with the density $`\rho _t`$ of Eq. (60), we sample first in $`\varphi `$ with a density proportional to $`\mathrm{ln}^2(\pi \gamma _{2t}/|\varphi |)`$. Then, with $`\varphi `$ chosen, we sample in $`A_{}`$ with a density proportional to $`\mathrm{log}(\gamma _{2t}/\omega )/\omega `$. Finally, having chosen $`\varphi `$ and $`A_{}`$, we sample in $`A_+`$ with a density proportional to the right hand side of Eq. (60). Taking into account the normalization conditions, this gives the result in Eq. (60) with

$`{\displaystyle \frac{1}{𝒜}}`$ $`=`$ $`{\displaystyle \frac{4\pi }{S_+\beta }}[\mathrm{ln}^2(\gamma _{2t})+2\mathrm{ln}(\gamma _{2t})+2]{\displaystyle \frac{1}{\gamma _{2t}\omega }}`$ (63)
$`\times \left\{1+{\displaystyle \frac{1}{2\mathrm{ln}(\gamma _{2t}/\omega )}}\mathrm{ln}\left({\displaystyle \frac{S_+1+\beta _{2t}S_+\omega ^2}{S_+1+\beta _{2t}\gamma _{2t}S_+\omega }}\right)\right\}`$
$`\times \left\{1{\displaystyle \frac{1}{2\mathrm{ln}^2\left(\pi \gamma _{2t}/|\varphi |\right)}}\left[\mathrm{ln}^2\left({\displaystyle \frac{\gamma _{2t}}{1+S_{}+|\varphi |/\pi }}\right)+\mathrm{ln}^2\left({\displaystyle \frac{\gamma _{2t}}{1S_{}+|\varphi |/\pi }}\right)\right]\right\}.`$
The function $`𝒜`$ may be a bit complicated, but it is quite benign. In particular, $`𝒜`$ is not singular or zero anywhere provided that $`\gamma _{2t}>\omega `$ holds everywhere. This requires that $`\gamma _{2t}>3`$.

## VIII Sampling for 2 to 3 scattering

In this section, we consider the sampling method for the third type of scattering singularity surface. We need to choose points $`\stackrel{}{l}\stackrel{}{l}_1`$ appropriate to the following case:

> A virtual parton with momentum $`\stackrel{}{l}`$ collides with a virtual parton with momentum $`\stackrel{}{l}`$ to produce a final state with partons carrying momenta $`\stackrel{}{l}_2`$, $`\stackrel{}{l}_3`$, and $`\stackrel{}{l}_2\stackrel{}{l}_3`$.

In this case, the scattering surface is the sphere $`2|\stackrel{}{l}|=|\stackrel{}{l}_2|+|\stackrel{}{l}_3|+|\stackrel{}{l}_2+\stackrel{}{l}_3|`$. The density of points should be large on this surface. There is no special point embedded in the surface where the density of points needs to be singular. However, the density should be enhanced near the point of this sphere where it intersects the ellipsoid in Fig. 6 or the narrower of the two ellipsoids in Fig. 8.

Since the scattering singularity surface is a sphere, we use spherical polar coordinates $`\{r,\mathrm{cos}\theta ,\varphi \}`$ to describe $`\stackrel{}{l}`$. We choose the $`\theta =0`$ axis along a certain direction $`\stackrel{}{P}_C`$. In the case (d) illustrated in Fig. 5, we define $`\stackrel{}{P}_C=\stackrel{}{p}_3`$. In the case (e) illustrated in Fig. 7, we define $`\stackrel{}{P}_C=\stackrel{}{p}_1`$ if $`|\stackrel{}{p}_1|>|\stackrel{}{p}_3|`$ or $`\stackrel{}{P}_C=\stackrel{}{p}_3`$ if $`|\stackrel{}{p}_3|>|\stackrel{}{p}_1|`$.

The jacobian of the transformation from $`\stackrel{}{l}`$ to $`\{r,\mathrm{cos}\theta ,\varphi \}`$ is

$$d\stackrel{}{l}=\frac{drd\mathrm{cos}\theta d\varphi }{\rho _{r\theta \varphi }},$$
(64)
where

$$\rho _{r\theta \varphi }=\frac{1}{r^2}.$$
(65)
If we sample points in the variables $`\{r,\mathrm{cos}\theta ,\varphi \}`$ with a density

$$\rho ^{}=\frac{dN}{drd\mathrm{cos}\theta d\varphi },$$
(66)
then the total density of points will be

$$\rho =\rho ^{}\times \rho _{r\theta \varphi }.$$
(67)
We now address the question of how to choose $`\rho ^{}`$. Our analysis follows closely that in Sec. VII B. The main idea is that there is a scattering singularity surface at $`r=S`$, where

$$S=\frac{1}{2}(|\stackrel{}{l}_2|+|\stackrel{}{l}_3|+|\stackrel{}{l}_2+\stackrel{}{l}_3|).$$
(68)
Thus the integrand has a factor

$$\frac{1}{(r/S)1+i\eta },$$
(69)
where $`\eta `$ is the amount of deformation of the $`r/S`$ contour. The amount of contour deformation vanishes quadratically as $`\{r,\mathrm{cos}\theta \}`$ approaches $`\{S_C,1\}`$, where

$$S_C=|\stackrel{}{P}_C|.$$
(70)
The point $`\{S_C,1\}`$ is the point $`\stackrel{}{l}=\stackrel{}{P}_C`$. That is, $`\stackrel{}{l}=\stackrel{}{p}_3`$ in Fig. 6 or in Fig. 8 when $`|\stackrel{}{p}_3|>|\stackrel{}{p}_1|`$. Thus, on the surface $`r=S`$, the amount of deformation $`\eta `$ can be estimated as $`\mathrm{\Delta }^2`$ where

$$\mathrm{\Delta }=\frac{1}{3}\sqrt{(1S_C/S)^2+(1\mathrm{cos}\theta )}.$$
(71)
Thus the absolute value of the integrand has a factor that can be estimated as

$$\frac{1}{|(r/S)1|}$$
(72)
for $`\mathrm{\Delta }1`$ and $`\mathrm{\Delta }^2|(r/S)1|\mathrm{\Delta }`$. We want $`\rho ^{}`$ to have a singularity of the same nature, so that the integrand divided by the density of points is singularity free. Furthermore, we would like $`\rho ^{}`$ to have an extra factor of $`1/\mathrm{\Delta }^2`$ to cancel the factor from the singularity associated with the narrow ellipsoid in Fig. 5 as we approach $`\theta =0`$. Thus we take

$$\rho ^{}=\frac{}{\mathrm{\Delta }\left(|(r/S)1|+\alpha _3\mathrm{\Delta }^2\right)\left(|(r/S)1|+\alpha _3\mathrm{\Delta }\right)},$$
(73)
where $`\alpha _3`$ is a fixed parameter and where $``$ is a rather complicated but nonsingular function that gets the normalization right:

$$=\frac{\alpha _3(1\mathrm{\Delta })}{18\pi S}\left[\mathrm{ln}\left(\frac{1}{\mathrm{\Delta }^2}\frac{1+\alpha _3\mathrm{\Delta }^2}{1+\alpha _3\mathrm{\Delta }}\right)\right]^1\left[\mathrm{ln}\left(1+\frac{2}{(1S_C/S)^2}\right)\right]^1.$$
(74)
## IX Sampling for 2 to 1 scattering

In this section, we consider the sampling method for the fourth type of scattering singularity surface:

> A virtual parton with momentum $`\stackrel{}{l}`$ combines with a virtual parton with momentum $`\stackrel{}{l}_2\stackrel{}{l}`$ to produce an on shell parton with momentum $`\stackrel{}{l}_2`$ that enters the final state.

As mentioned in Sec. IV, the 2 to 1 scattering singularity surface is exceptional in that there is no singularity except for the two singular points at $`\stackrel{}{l}=0`$ and $`\stackrel{}{l}_2\stackrel{}{l}=0`$. Typically a 2 to 1 scattering subdiagram is part of a 2 to 2 (s) or 2 to 2 (t) scattering subdiagram and the two singularities are provided for by the 2 to 2 (s) or 2 to 2 (t) sampling methods. The exception is in the case of a self-energy subgraph that is connected to a final state parton. In this case, the 2 to 2 (t), 2 to 2 (s), and 2 to 3 sampling methods do not apply and we need an explicit 2 to 1 sampling method. Thus we apply a 2 to 1 sampling method only in the case of a self-energy subgraph connected to a final state parton.

We will choose $`\stackrel{}{l}`$ with a density that is a linear combination of five densities:

$$\rho (\stackrel{}{l})=\alpha _{1a}\rho _a(\stackrel{}{l})+\alpha _{1b}\rho _b(\stackrel{}{l})+\alpha _{1a}\rho _a(\stackrel{}{l}\stackrel{}{l}_2)+\alpha _{1b}\rho _b(\stackrel{}{l}\stackrel{}{l}_2)+(12\alpha _{1a}2\alpha _{1b})\rho _c(\stackrel{}{l}).$$
(75)
Here $`\alpha _{1a}`$ and $`\alpha _{1b}`$ are fixed positive parameters with $`(12\alpha _{1a}2\alpha _{1b})`$ also positive.

For $`\rho _a(\stackrel{}{l})`$ we take

$$\rho _a(\stackrel{}{l})=\frac{|\stackrel{}{l}_2|}{4\pi }\frac{1}{\stackrel{}{l}^{\mathrm{\hspace{0.17em}2}}}\frac{1}{\left(|\stackrel{}{l}|+|\stackrel{}{l}_2|\right)^2}.$$
(76)
The density $`\rho _b`$ is a simple variation on this with $`|\stackrel{}{l}_2|`$ replaced by a scale $`M`$:

$$\rho _b(\stackrel{}{l})=\frac{M}{4\pi }\frac{1}{\stackrel{}{l}^{\mathrm{\hspace{0.17em}2}}}\frac{1}{\left(|\stackrel{}{l}|+M\right)^2}.$$
(77)
Here $`M=(|\stackrel{}{l}_2|+|\stackrel{}{l}_3|+|\stackrel{}{l}_2+\stackrel{}{l}_3|)/3`$, where the momenta of the particles in the final state are $`\stackrel{}{l}_2`$, $`\stackrel{}{l}_3`$, and $`\stackrel{}{l}_2\stackrel{}{l}_3`$. The $`\rho _a`$ and $`\rho _b`$ terms provide singularities at $`\stackrel{}{l}^2=0`$ and $`(\stackrel{}{l}_2\stackrel{}{l})^2=0`$.

The $`\rho _c`$ term provides a non-singular density of points in the neighborhood of the collinear line $`\stackrel{}{l}=x\stackrel{}{l}_2`$ with $`0<x<1`$. For $`\rho _c(\stackrel{}{l})`$ we take

$$\rho _c(\stackrel{}{l})=\frac{\beta _1\gamma _1^2|\stackrel{}{l}_2|}{2\pi }\frac{1}{\sqrt{\beta _1^2+(x1/2)^2}\left(\sqrt{\beta _1^2+(x1/2)^2}+\beta _1\right)\left(\stackrel{}{l}_T^{\mathrm{\hspace{0.17em}2}}+\gamma _1^2\stackrel{}{l}_2^{\mathrm{\hspace{0.17em}2}}\right)^2}.$$
(78)
Here $`\beta _1`$ and $`\gamma _1`$ are fixed parameters and $`x`$ and $`\stackrel{}{l}_T`$ are defined by writing

$$\stackrel{}{l}=x\stackrel{}{l}_2+\stackrel{}{l}_T,$$
(79)
where $`\stackrel{}{l}_T`$ is the part of $`\stackrel{}{l}`$ that is orthogonal to $`\stackrel{}{l}_2`$.

## X Conclusions

In a next-to-leading order calculation of three jet cross sections and similar observables in electron-positron annihilation, a given graph leads to several contributions to the cross section. Each of these contributions has the form of a measurement function times a quantum amplitude times a complex conjugate amplitude, as in Fig. 1. We must integrate the sum of these contributions over three loop momenta. In this paper, I have presented a method for sampling the space of loop momenta in order to perform the integrations by the Monte Carlo method.

The main organizing principle for this sampling is to construct the density of integration points as a sum and to have one or more terms in this sum correspond to each possible three parton final state for the graph.

A cut graph with a three parton final state has a virtual loop. We are thus led to consider a simple problem in quantum mechanics: where in the space of the loop momentum are the singularities for two body to $`n`$ body scatterings? The generic answer is that the singularities lie on ellipsoidal surfaces. This leads us to choose elliptical/hyperbolic coordinates for the loop momentum space. For the most part, we avoid letting the integrand be truly singular at these surfaces by giving the loop momentum a suitable imaginary part according to the contour deformation recipe of Ref. . Nevertheless, it is helpful to make the density of integration points large near these ellipsoidal surfaces.

We are left with special points on the scattering singularity surfaces where the contour deformation vanishes, so that the integrand is actually singular. This singularity corresponds to soft parton exchange and is significant in gauge field theories. We have seen how to give the density of points a singularity structure that matches that of the integrand as one approaches a soft singularity.

The sampling method that has been described here is far from optimal. It is, however, at least reasonably systematic, and it gives good results for three jet observables that do not get significant contributions from parton final states that are near to two jet configurations.

###### Acknowledgements.

This work was supported in part by the U.S. Department of Energy.


text:
# Comment on “Lattice QCD analysis of the strangeness magnetic moment of the nucleon”

\[

## Abstract

The recent chirally extrapolated result of Leinweber and Thomas \[Phys. Rev. D 62, 074505 (2000), or LT\] for the nucleon strangeness form factor $`G_M^s(0)`$ = -0.16 $`\pm `$ 0.18 $`\mu _N`$ differs markedly from the earlier result -0.75 $`\pm `$ 0.30 $`\mu _N`$ obtained by Leinweber \[Phys. Rev. D 53, 5115 (1996)\] from the same lattice data. An unresolved problem in the LT analysis of lattice data is identified and addressed. A value of $`G_M^s(0)`$ = -0.55 $`\pm `$ 0.37 $`\mu _N`$ is obtained at $`R_d^s`$ = 0.55 by extrapolating only the nucleon isoscalar lattice data.

\]

Past lattice calculations made by Leinweber and collaborators  of baryon magnetic moments (MMs) have been restricted to massive quarks equivalent to pion masses $`m`$ greater than 600 MeV. The calculated moments have to be extrapolated to the physical pion mass $`m_\pi `$ = 140 MeV before comparison with experiment. Extrapolating linearly in $`m^2`$, Leinweber  has estimated the nucleon strangeness form factor at momentum transfer $`Q^2=0`$ to be $`G_M^s(0)`$ = -0.75(30) in units of $`\mu _N`$. (The number inside parentheses gives the statistical uncertainty in the last digits.) A more sophisticated nonlinear extrapolation of the same lattice data using a chiral extrapolation ($`\chi `$E) formula proposed by Leinweber, Lu and Thomas (LLT)  has recently been made by Leinweber and Thomas (LT) . The result is -0.16(18).

This LT estimate contains an unresolved problem caused by the use of an extrapolated valence $`u`$-quark contribution $`u_{\mathrm{\Xi }^0}`$ for the baryon $`\mathrm{\Xi }^0`$ that disagrees with experiment to renormalize another extrapolated term $`u_n`$. The purpose of this Comment is to suggest that the lattice data should be renormalized before the $`\chi `$E, and to point out that $`G_M^s(0)`$ can be extracted more readily from the nucleon isoscalar lattice data. The result is -0.55(37), where the error comes from uncertainties in the extrapolation formula, in the least-square fitting, and in the SU(6) constants $`F`$ and $`D`$.

To extrapolate a MM term $`\mu (m)`$ calculated on a lattice for several large values of the pion mass $`m`$ to its physical value at $`m_\pi `$, LT use a Padé approximant given by their Eq. (13) (denoted here LT13) that is the inverse of a quadratic Taylor approximant in $`m`$. The extrapolated results, needed in this Comment, are read from Figs. 4 and 5 of  and given in Table I. They agree well with experimental values  where known, as pointed out in  for $`\mu _p=p`$ and $`\mu _n=n`$.

For technical reasons, LT use a simple quadratic Taylor approximant in $`m`$, Eq. (14) in LT and called here LT14, for $`u_{\mathrm{\Xi }^0}`$. (In terms of $`s`$ = $`m^2`$, LT14 is the often-used linear approximant improved by a certain leading nonanalytic (LNA) term $`\chi m`$.) Table I shows that the extrapolated value for $`u_{\mathrm{\Xi }^0}`$ differs significantly from the experimental value $`\mathrm{\Xi }^0\mathrm{\Xi }^{}`$. A concise review of these results in greater details has been given by .

The discrepancy from experiment is not due to the use of LT14, as I have explicitly obtained an LT13 fit giving the same extrapolated value of $`u_{\mathrm{\Xi }^0}`$ = -0.37 (and -0.38 for LT14). The discrepancy comes instead from the sign of $`\chi `$ that controls the curvature direction caused by the LNA term. In the LT model, $`\chi `$ has opposite signs for $`u_n`$ and $`u_{\mathrm{\Xi }^0}`$, because the dominant pion cloud has opposite signs, namely $`n`$ $`p\pi ^{}n`$ compared to $`\mathrm{\Xi }^0`$ $``$ $`\mathrm{\Xi }^{}\pi ^+`$ $`\mathrm{\Xi }^0`$ . This sign makes it difficult for the $`\chi `$E of $`u_{\mathrm{\Xi }^0}`$ to curve towards the experimental point. This is the unresolved problem in LT mentioned at the beginning of this Comment

To extract $`G_M^s`$ (dropping the $`Q^2`$ argument for notational simplicity), LT use a certain “ratio” method of Leinweber  under which the disconnected-loop (DL) contribution $`\mu _N^{\mathrm{DL}}`$ of the nucleon sea is calculated from one of the following formulas:

$`\mu _N^{\mathrm{DL}}={\displaystyle \frac{1}{3}}(p+2nf_nu_n),f_n={\displaystyle \frac{\mathrm{\Xi }^0\mathrm{\Xi }^{}}{u_{\mathrm{\Xi }^0}}};`$ (1)
$`\mu _N^{\mathrm{DL}}={\displaystyle \frac{1}{3}}(2p+nf_pu_p),f_p={\displaystyle \frac{\mathrm{\Sigma }^+\mathrm{\Sigma }^{}}{u_{\mathrm{\Sigma }^+}}}.`$ (2)
That is, the extrapolated value of $`u_n`$ is not used directly, but only after renormalization by $`f_n`$. Whenever $`u_{\mathrm{\Xi }^0}`$ extrapolates correctly to its experimental value $`\mathrm{\Xi }^0\mathrm{\Xi }^{}`$, the correction factor $`f_n`$ is 1 and therefore has no effect. The correction factor differs from 1 only when the extrapolation is unsuccessful. Then the argument for using the nontrivial correction factor $`f_n1`$ is that the ratio $`u_n/u_{\mathrm{\Xi }^0}`$ that appears has smaller systematic errors .

The large renormalization $`f_n`$ = 1.62 found for Eq. (1) is the direct consequence of the aformentioned failure of the LT extrapolation for $`u_{\mathrm{\Xi }^0}`$. In contrast, $`f_p`$ = 0.973 needed in Eq. (2) contains little renormalization because LT13 works well there. The LT renormalized values (actually $`f_iu_i`$) are shown in Table I in the column marked LTR with errors from LT.

Is this postulated cancellation of systematic errors after LT extrapolations real? On the lattice, $`u_n`$ and $`u_{\mathrm{\Xi }^0}`$ have very similar values, and their ratios are likely to contain smaller systematic errors. Unfortunately, these ratios cannot be used directly in the $`\chi `$Es, and the $`u`$’s have to be extrapolated separately for the lattice data used. Hence the correct procedure is to renormalize one set of lattice data to fit an available experimental MM (thus ensuring that $`f_i`$ = 1 always), and to use the same renormalization on the second set of lattice data before the $`\chi `$E.

This scheme can readily be executed for $`u_{\mathrm{\Sigma }^+}`$ since the original extrapolation is already close to experiment. With the lattice data for different $`m`$ actually correlated, moving up and down together, one can simply take $`\{u_i+x\mathrm{\Delta }u_i\}`$ as the corrected input with the unknown $`x`$ chosen to be -0.23 to make the extrapolated value of $`u_{\mathrm{\Sigma }^+}`$ fit $`\mathrm{\Sigma }^+\mathrm{\Sigma }^{}`$. The same $`x`$ is then used to renormalize the lattice data for $`u_p`$. The LT13 extrapolated result, $`u_p`$ = 4.16, agrees with that from Eq. (2).

When the same procedure is applied to $`u_{\mathrm{\Xi }^0}`$, it reaches a minimum of only -0.52, a little short of the experimental value. (The solution for the input $`y\{u_i\}`$ scaled by an adjustable factor $`y`$ reaches down to only -0.47.) There are LT14 solutions ($`u_p`$ = 2.34, $`u_n`$ = -1.54) that are very different from those of the Leinweber equations. Since LT14 also does very poorly for the nucleon MMs themselves, giving $`p`$ = 3.14 and $`n`$ = -2.39, it appears to be unreliable.

The lattice data can also be renormalized by tilting and bending. This is conveniently realized by including the experimental point in an LT13 fit and using the resulting set of ratios $`\{u_i(\mathrm{output})/u_i(\mathrm{input})\}`$ to renormalize the second set of lattice data before its $`\chi `$E. The resulting renormalized lattice data (RLD) extrapolate to $`u_n`$ = -0.64(15), showing only a quarter of the renormalization effect from Eq. (1), while $`u_p`$ = 4.13(24) agrees well with Eq. (2). These results suggest that Eq. (2) is reliable, but Eq. (1) greatly over-estimates the renormalization. The statistical errors shown are obtained by treating the lattice data as correlated, i.e., calculated from $`\{u_i\pm p\mathrm{\Delta }u_i\}`$, where $`p`$ = $`\sqrt{1(N_p/N_d)}`$ = $`1/\sqrt{3}`$, where $`N_d`$ = 3 is the number of lattice data and $`N_p`$ = 2 is the number of fitted parameters .

An alternative explanation of the discrepancy in $`u_{\mathrm{\Xi }^0}`$ that cannot be excluded at the present time is that its $`\chi `$ parameter is different from that calculated in LT. For example, the $`\chi `$ parameter needed to fit the experimental point together with the original lattice data is $`\chi (u_{\mathrm{\Xi }^0})`$ = 0.2 (-0.3/+0.5), showing the change of sign needed for the $`\chi `$E to curve towards the experimental point. In other words, a relatively small change in the LT value of $`\chi (u_{\mathrm{\Xi }^0})`$ = -0.4 is enough to change its sign. In the rest of this Comment, I shall work only with the LT model.

The value of $`G_M^s`$ = -0.16(18) reported by LT is obtained from Eq.(1) using the ratio $`u_n/u_{\mathrm{\Xi }^0}`$ = 1.51(37). (The reported statistical error of 0.18 for $`G_M^s`$ is a misprint. It should read 0.27.) LT has also obtained $`G_M^s`$ = -0.57(42) from Eq.(2) using the ratio $`u_p/u_{\mathrm{\Sigma }^+}`$ = 1.14(11), unchanged from the Leinweber linearly extrapolated value of 1.14(8). This is the better estimate because the renormalization involved is so much smaller and better established.

The second issue raised in this Comment is that the nucleon isoscalar quantity $`G_M^s`$ can be extracted from nucleon isoscalar lattice data alone without using any information from $`u_{\mathrm{\Xi }^0}`$ or $`u_{\mathrm{\Sigma }^+}`$. The resulting physical picture for $`G_M^s`$ is sufficiently simple to permit a rather unique answer to be obtained for any extrapolation formula used such as LT13. Rough estimates of other uncertainties can then be made readily.

The quantity to be extracted is the isoscalar DL contribution that, under isospin symmetry, is simply

$`\mu _N^{\mathrm{DL}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(p+n){\displaystyle \frac{1}{6}}(u_p+u_n).`$ (3)
The numerical answer is $`\mu _N^{\mathrm{DL}}`$ = -0.18(5) using the experimental value of $`p+n`$ and the unrenormalized LT results shown in Table I. The answer is -0.10(7) with the original LT renormalization, and -0.14(9) using RLD to renormalize $`u_p`$ and $`u_n`$. The statistical errors shown are obtained from those in $`u_p`$ and $`u_n`$ treated as independent.

To do better in both value and error, it is necessary to extrapolate $`u_p`$ \+ $`u_n`$ together. Fortunately, there is enough information in the lattice data tabulated in LWD  to extract the covariance $`(\mathrm{\Delta }u_p)(\mathrm{\Delta }u_n)`$ and to calculate the statistical error for any linear combination $`au_1`$ \+ $`bu_2`$ on the lattice. The combination is then extrapolated to $`m_\pi `$ by using the chiral parameter $`a\chi _1`$ \+ $`b\chi _2`$.

The least-square fitting program and input lattice data used here are first validated by checking against published results. The top panel of Table I gives a comparison with the results of LLT  for the nucleon MMs $`p`$ and $`n`$. The data used come from two independent sources . My errors are obtained by assuming that all lattice data are correlated (i.e., using $`N_p`$ = 2 and $`N_d`$ = 6). The chiral parameters $`\chi `$ are defined in terms of the SU(6) constants $`F`$ and $`D`$ . The one-loop corrected values  are used in both LLT and LT. The table shows good agreement in both extrapolated values and fitting errors.

In the second panel of the table, extrapolations are obtained for the correlated lattice data  used in LT. Both values and errors agree with their results.

The third panel of Table I shows the best results from Tabel II obtained by extrapolating $`u_p`$ $`\pm `$ $`u_n`$ together, using different $`\chi `$ parameters for flavor and valence contributions. The quality of this single-step extrapolation is first checked for the nucleon MMs. The isovector moment $`pn`$ = $`u_p`$ \- $`u_n`$ on the lattice because the isoscalar sea contributes nothing. Table II shows that the $`pn`$ result of 4.65 for LT13 agrees well with the value of 4.68 from separate extrapolations. Both differ somewhat from the value $`u_p`$ \- $`u_n`$ = 4.81 from separate extrapolations. The isoscalar moment extrapolates to $`p+n`$ = 1.03, in agreement with the value of 1.10 from separate extrapolations and with experiment. In contrast, LT14 appears inadequate for both nucleon MMs.

The much more successful formula LT13 proposed by LLT has the correct $`1/m_q1/m^2`$ behavior for heavy quarks. Its comparative success in chiral extrapolations comes from a reduction of the effect of the LNA term when placed in the denominator. Finally, the importance of having the correct heavy-quark limit is illustrated by using another inverse Taylor approximant M4D where the quadratic term in $`m`$ in the denominator is replaced by an incorrect quartic term. It works well too for the large isovector moment $`(pn)/2`$ that dominates the nucleon MMs. This suggests once again that the LNA term should be in the denominator.

The extrapolated value of $`u_p+u_n`$ is also needed to calculate $`\mu _N^{\mathrm{DL}}`$ from Eq. (3). It is obtained from the lattice data for $`u_p`$ \+ $`u_n`$ by extrapolating with a valence parameter $`\chi ^v`$ different from the flavor value $`\chi ^f`$ = 0 used for $`p`$ \+ $`n`$. The results obtained with different extrapolation formulas are shown in Table II. The very small statistical error found for $`\mu _N^{\mathrm{DL}}`$, denoted in the table as $`\mu _N^{\mathrm{DL}}(1)`$, comes from the fact that while each term in Eq. (3) varies substantially as all lattice data move up and down together, the difference between the two terms remains essentially unchanged.

To determine if $`p+n`$ should be fitted in extracting $`\mu _N^{\mathrm{DL}}`$, I renormalize the lattice input data to $`\{u_i+x\mathrm{\Delta }u_i\}`$ where $`x`$ is adjusted to reproduce the experimental value. The changes for all three formulas, also given in the table, are all quite small, showing that the procedure is quite robust.

Unfortunately, there is a problem in the extrapolations for $`p`$ \+ $`n`$, and also for $`p`$ and $`n`$ separately. The lattice data used do not contain any DL contribution, which has in fact been assumed negligibly small in all LT extrapolations. Unpublished lattice data from Dong, Liu and Williams(DLW)  seem to indicate that the DL contribution is not necessarily negligible in the mass range $`m`$ 0.6-0.9 GeV of the LWD data. Unfortunately, the covariances for these DL lattice data have not been saved, thus precluding a quantitative estimate of the missing DL contribution for different $`m`$’s. This missing DL term is formally twice as strong and fractionally 4-6 times more important in $`p`$ \+ $`n`$ as it is in each MM separately. It gives an unknown systematic error to the extrapolated $`p`$ \+ $`n`$.

Given this unresolved problem in the lattice data, one might want to use the experimental value of $`p`$ \+ $`n`$ in Eq. (3). The results, shown in the table as $`\mu _N^{\mathrm{DL}}(2)`$, can differ significantly from $`\mu _N^{\mathrm{DL}}(1)`$. The large error shown comes from the statistical fitting error of $`u_p`$ \+ $`u_n`$.

These two estimates of $`\mu _N^{\mathrm{DL}}`$ agree if the extrapolated $`u_p`$ \+ $`u_n`$ agrees with experiment. It is therefore interesting to see what happens when the lattice data for $`u_{\mathrm{\Xi }^0}`$ and $`u_{\mathrm{\Sigma }^+}`$ are used to renormalize the lattice data for $`u_n`$ and $`u_p`$ before extrapolating $`u_p`$ $`\pm `$ $`u_n`$ by LT13 and M4D. The results, given in the lower panel of Table II, show that the extrapolated $`p`$ \+ $`n`$ by LT13 is now close to the experimental value. I take $`\mu _N^{\mathrm{DL}}(2)`$ as the better estimate because it contains no systematic error from $`p`$ \+ $`n`$ and a more conservative error estimate. This is the “best” result shown in Table I. The difference of 0.07 between LT13 and M4D shall be taken to be an additional uncertainty arising from incomplete knowledge of the correct $`\chi `$E formula. This large uncertainty has not been included in LT.

Errors also arise from the choice of the SU(6) constants $`F`$ and $`D`$ on which the $`\chi `$s depend. The $`\chi `$ parameter for the isoscalar MM $`p+n`$ is always $`\chi ^f`$ = 0 for any $`F`$ and $`D`$. The valence parameter $`\chi ^v`$ for $`u_p+u_n`$ is proportional to (and has the same sign as)

$`\beta ^v=a(F+D)^2,a={\displaystyle \frac{2}{3}}\left[{\displaystyle \frac{56r+9r^2}{(1+r)^2}}\right],`$ (4)
where $`rF/D`$. The factor $`a`$ is always positive, so that both $`\mu _N^{\mathrm{DL}}`$ and $`G_M^s`$ extracted from it are always negative in LT extrapolations. The factor $`a`$ has a minimum of 2.39 at $`r`$ = 0.74, and the values 2.40 and 2.45, respectively, at the one-loop corrected value of $`r`$ = 0.66 and the tree-level value of $`r`$ = 0.58 . Hence the error from uncertainties in $`r`$ is small and will be neglected.

With $`r`$ = 0.66, the uncertainty in $`F+D`$ can be determined by treating it as a free parameter in LT13 fitted to both experimental and lattice data for each nucleon MM $`p`$ or $`n`$. The expanded data set used by LLT then gives

$`(F+D)_p=0.96(20),(F+D)_n=1.00(20),`$ (5)
$`\mathrm{\Delta }(u_p+u_n)=0.21,\mathrm{\Delta }\mu _N^{\mathrm{DL}}=0.04.`$ (6)
This yields a total error of 0.10 for $`\mu _N^{\mathrm{DL}}`$ reported in Table I when quadratically combined with the errors of 0.06 (from least-square fitting) and 0.07 (from the $`\chi `$E formula).

The resulting values of $`G_M^s`$ are shown in the lowest panel of Table I for $`R_d^s`$ = 0.55 and 0.59 as used by LT and LLT, respectively. The value of -0.55(37) at $`R_d^s`$ = 0.55 agrees with the result -0.36(20) of DLW  obtained from a direct evaluation of the DL contribution. The value $`R_d^s`$ = 0.55 used here is actually their ratio of directly evaluated strange to $`u/d`$ DL contributions.

From the perspective of Eq. (3), the Leinweber Eq. (1) or (2) contains a spurious nucleon isovector term

$`\mathrm{\Delta }\mu _N^{\mathrm{DL}}={\displaystyle \frac{1}{6}}[pnf_i(u_pu_n)],`$ (7)
where $`f_i`$ = $`f_p`$ or $`f_n`$. With $`f_p`$ close to 1, this term is very small in Eq. (2). This is another reason why Eq. (2) is better than Eq. (1).

Additional systematic errors not taken into consideration here include effects from finite-volume and finite-spacing errors and from the quenched approximation not already simulated by the $`\chi `$E formula used . They also include uncertainties in the ratio $`R_d^s`$ and differences between the LT model and true QCD, such as the contributions of neglected kaon loops and of $`m_q\mathrm{ln}m_q`$ terms, in as far as they affect the parameter analogous to $`\beta ^v`$ of Eq. (4) that drives the DL contribution in the $`\chi `$E.

I would like to thank Drs. D. Leinweber and K.F. Liu for helpful comments.