X Artru
email: x.artru@ipnl.in2p3.fr

C Ray
email: c.ray@ipnl.in2p3.fr

Electron acceleration in a helical waveguide 1

Travelling wave in a helical wave guide is considered for electron acceleration. A first determination of the travelling wave modes using a partial wave expansion (PWE) and a point matching method is presented. It gives a rapid solution for moderate deformation of the guide relative to a straight cylinder. Strong deformations will give higher accelerating gradient but the PWE diverges. Two methods overcoming this difficulties are suggested.

Introduction

An electromagnetic wave in a straight wave guide cannot accelerate particles because its phase velocity v ph is larger than c. In conventional linear accelerators, v ph is reduced by a chain of resonant cavities along the guide. A scheme not using a resonance phenomenon is the travelling wave tube (TWT). One type of TWT is a cylindrical tube with a conducting helical wire inside [START_REF] Klein | Proceedings of the 1972 proton linear accelerator conference[END_REF][START_REF] Aron | Stanford[END_REF]. Another type is a tube with a thick dielectric wall [START_REF] Zhang | [END_REF].

In this paper we study another kind of helical device : the helically bent tube. The main "photon path" is the helix passing trough the centres of the local cross sections of the tube. For a strong enough bending, the difference in length between this path and a straight path is such that the effective v ph is smaller than c, although the local v ph is larger than c.

One may expect several advantages in using helically bent tube in place of accelerating cavities : easy construction, easy cooling, easy maintainance of the wall surface using an helical arm, large frequency spacing of the modes, large group velocity which allow a fast filling of the guide by the wave. Alternatively, the wave can be generated by a driving beam. Like with the dielectric wall in [4], several modes can cooperate to transfer the energy from one beam to the other.

The helical wave guide is defined by its axis Oz, its period Λ = 2π/q and the base curve, which is the section of the wall with the z = 0 plane, of parametric equation r b (s) = (x b (s), y b (s)). The point (x, y) = (0, 0) is supposed to be inside the base. The wall surface is generated by helices H s of parametric equation

(x, y, z) = (R[qz]r b (s), z) (1)
where the 2×2 matrix R[qz] describes the rotation of angle qz about 0z. One method to make a helical waveguide is to draw helical grooves on the internal face of an initially cylindral tube. The case of rectangular grooves has been studied in [5,[START_REF] Wallet | Technical characteristics of a novel helical-groove traveling-wave tube structure[END_REF]. The corresponding wall base curves have indentations. Here we will restrict to base curves with smooth curvatures. This work is the continuation of a first study [START_REF] Artru | Proceedings of Channeling 2008[END_REF], where the electromagnetic wave modes were calculated using a partial wave expansion (PWE) associated to a point matching (PM) method. Since that time, we realized that PWE is limited to weak enough helical bending and sufficiently "round" base curves ; if it is not the case, the PWE converges only inside a cylinder which does not contain the whole guide.

In this paper we will recall the PWE-PM method, some of our old results, present some new ones and discuss the limitations of this method. Two alternative methods will be evoked.

Propagating modes in the scalar field case

We first present the simpler case of a scalar field Φ(t, x, y, z), which obeys the Klein-Gordon equation (in units c = 1)

(-∂ 2 t + ∂ 2 x + ∂ 2 y + ∂ 2 z ) Φ(t, x, y, z) = 0 , (2)
with the boundary condition Φ(t, x, y, z) = 0 on the wall. These two equations are invariant under a time translation and under a helical displacement {translation by ∆z} × {rotation by q∆z} .

In the transverse plane we introduce a "rotating frame", the coordinates (X, Y ) of which are related to the fixed frame ones by x + iy = e iqz (X + iY ). On can therefore look for a solution of the form

Φ(t, x, y, z) = e iP z-iωt Ψ(X, Y ) (4)
which has definite reduced frequencies ω and definite eigenvalue P of the operator

P helic = p z + qJ z = -i∂/∂ z -iq∂/∂ ϕ . (5)
2.1 The partial wave expansion (scalar field case)

Let us temporarily assume that a mode can be expanded in free cylindrical partial waves which are eigenstates of J z and p z with eigenvalues l and p respectively. A single partial wave writes Φ ω,p,l (t, x, y, z) = e -iωt e ipz e ilϕ ψ l,k (r) ,

where (r, ϕ, z) are the cylindrical coordinates and

ψ l,k (r) = J |l| (kr)/k |l| = I |l| (κr)/κ |l| . (7)
J l and I l are the regular Bessel functions of first and second kind, k = (ω 2 -p 2 + i0) 1/2 is the radial momentum and κ ≡ -ik = (p 2 -ω 2 -i0) a l e ip(l)z+ilϕ ψ l,k(l) (r) ,

with p(l) or p l = P -lq , k(l) or

k l = ω 2 -p 2 l + i0 1/2 , (9)
so that each partial wave is an eigenstate of P helic with eigenvalue P . It is sufficient to write the boundary condition at t = z = 0, thanks to the invariances under time translation and helical displacement. For every point r b (s), of polar coordinates (r b (s), ϕ b (s)), of the base curve one has the constraint

a l exp[ilϕ b (s)] ψ l,k [r b (s)] = 0 . ( 10
)
Figure 1: Matching points on the base curve.

The point matching method in the scalar field case

One keeps only a finite set L of partial waves and selects an equal number of matching points r 1 , ...r N of the base curve (see Fig. • the r n 's should give a precise enough definition of the base curve,

• L should contain l's of both signs, however with more l's of the same sign as q.

• l (ϕ n+1 -ϕ n ) for l ∈ L should not be too large compared to unity,

• k l (r n+1 -r n ) for l ∈ L should not be too large compared to unity.

The a l are approximatively determined by the system of linear equations

l∈L a l e ilϕn ψ l,k (r n ) = 0 , 1 ≤ n ≤ N , (11)
(r n , ϕ n ) being the polar coodinates of r n . A nonzero solution exists if the determinant est zero:

det{C nl } = 0 , where C nl = e ilϕn ψ l,k (r n ) . (12)
This equation, which relates ω and P , is the dispersion relation in the PWE-PM approximation. One may impose ω and look for the roots of (12) as a function of P , or vice-versa. Alternatively, in view of accelerating highly relativistic electrons one may impose ω/P ≡ v ph = v e ≃ 1 and look for the roots of (12) as a function of ω. Using a graphical software, the roots are easily seen as dips in the curve of ln |det{C nl }| versus the free variable. If (det{C nl }) is e iα times a real number, the roots are precisely located as jumps of the curve arg (ie -iα det{C nl }) from ±π/2 to ∓π/2. Provided that the partial wave expansion converges, the quality of the approximation is improved by increasing N. The convergence is obtained first for the lowest modes (which have smallest ω at fixed P or v ph , or largest P at fixed ω), then for higher and higher modes. A first criterium is the stability of the root when N is increased by one or two units, when L is translated by one or two units or when one moves the matching points. The ultimate criterium is the stability of the coefficients a l , at least for l not close to l inf or l sup .

Case of a symmetric base. If the base is symmetrical about the x-axis, one can choose the {r n } to be symmetrical by pairs, except for those which are on the x-axis. Then i ν det{C nl } is real, ν being the number of pairs. The a l can be taken all real, which gives Φ(0, x, -y, 0) = Φ * (0, x, y, 0) .

(

) 13
If the base is symmetrical about Oz, the modes are either even or odd under the parity symmetry (x, y) → (-x, -y), therefore L has to contain either even or odd l's. More generally, if the base has a rotational symmetry of order S, the l's should be equal modulo S. To avoid redundant equations in (12), the points r n should be all different modulo the symmetry. It suffices to choose them all on one sub-period of the base curve.

Numerical exemple with the scalar field. The simplest base curve is a circle or radius a which is de-centered by εa with respect to the z-axis. We took ε = 0.4, q = 1 and imposed the phase velocity v ph = 0.99999. The lowest modes are at ωa = 5.0 and 7.0. They are stable at the 0.1 precision when L contains at least the set [-1, +5]. Fig. 2 displays the cross section Ψ(X, Y ) of the field at z = t = 0, for ωa = 5.0 mode. Note the symmetry (13). The maximum of |Re Ψ| is not at the centre of the base, but farther from the helix axis. This is a manifestation of the centrifugal effect.

Propagating modes of the Maxwell field

We gather the fields E and B in a 6-component vector F = {E, B}. In analogy with (4) we look for solutions of the form

F(t, x, y, z) = e iP z-iωt R[qz] F (X, Y ) (14) with F (X, Y ) = {E(X, Y ), B(X, Y )}. R[qz]
is the rotation matrix of (1) applied to the vectors E and B. The physical field is the real part of F.

Boundary conditions. We introduce the vectors

σ(s) = r ′ b (s) =   x ′ b (s) y ′ b (s) 0   , τ (s) = ẑ + q ẑ × r b (s) =   -qy b (s) qx b (s) 1   , (15)
which are tangent to the wall at z = 0, and the normal vector ν(s) = σ(s) × τ (s). The boundary conditions are

   E • σ = 0 E • τ = 0 B • ν = 0 →    E r σ r + E ϕ σ ϕ = 0 qr E ϕ + E z = 0 B r σ ϕ + (qrB z -B ϕ )σ r = 0 (16)
The right set is written in cylindrical coordinates. These 3 conditions are not independent. Only two suffice.

Partial wave expansion. F can be expanded in free partial waves like Φ in [START_REF] Petit | [END_REF]. However, since the photon has two polarisation states, there are two waves for each l, which now denotes the total (spin + orbital) angular momentum. We take the states of right and left chirality, or photon helicity χ = ±1. The analogue of ψ l,k (r) in [START_REF] Wallet | Technical characteristics of a novel helical-groove traveling-wave tube structure[END_REF][START_REF] Petit | [END_REF] is F ω,p,l,χ (r) = E ω,p,l,χ (r), B ω,p,l,χ (r) . For l = 0, the cylindrical components are

E ω,p,0,χ (r) =   E r E ϕ E z   ω,p,0,χ =   -ip χ ψ 1 (r) ω ψ 1 (r) χ ψ 0 (r)   , B ω,p,0,χ (r) = -iχ E ω,p,0,χ (r) . (17) For l = 0, E ω,p,l,χ =   -i [ψ λ-1 + (p -χ ′ ω) 2 ψ λ+1 ] s l [ψ λ-1 -(p -χ ′ ω) 2 ψ λ+1 ] 2(p -χ ′ ω) ψ λ   , B ω,p,l,χ = -iχ E ω,p,l,χ . (18)
Here s l = sign(l), χ ′ = s l χ, λ = |l| and ψ λ = J λ (kr)/k λ . These states have finite and linearly independent limits when p 2 → ω 2 . The analogue of ( 8) is

F(t, x, y, z) = e -iωt l
e ilϕ e ip(l)z a l F ω,p(l),l,χ=+s l (r) + b l F ω,p(l),l,χ=-s l (r) .( 19)

We made the choice to associate a l and b l with positive and negative χ ′ instead of χ. For l = 0 we set s l = +1 in (19), but we can alternatively take the (TM, TE) basis as in [START_REF] Artru | Proceedings of Channeling 2008[END_REF].

The point matching method in the Maxwell field case. Generalising the PWE-PM method of Section 2.2 we keep only N values of l and write the first two conditions of (16) for N points r n of the base curve. It gives 2N linear equations linking N coefficients a l and N coefficients b l . The 2N × 2N matrix of the system is given in [START_REF] Artru | Proceedings of Channeling 2008[END_REF]. The dispersion relation is the vanishing of its determinant.

Numerical results for the Maxwell field. Taking the same de-centered circular base curve as in the scalar field case ε = 0.4, but q = 1 and imposing v ph = 0.999, we obtain the lowest modes at ωa = 2.28, 3.6, 4. Increasing the guide periods Λ = 2π/q while imposing fixed v ph < 1 makes the ω of each mode increase until the mode disappear at ω = ∞. For q = 0.7, the two lowest modes are at ωa = 2.7 and 4.2 with |E z (0)/B z (0)| = 0.25 and 0.26 respectively. The lowest mode has f.o.m ≃ 1/9.

With a stronger helical bending of the waveguide (larger ε and q), the mode frequencies decrease, but the convergence of the PWE-PM method becomes questionable. For ε = 0.6, q = 1.5 and v ph = 0.999, three first modes are seen at ωa = 2.2, 3.3 and 4.0 with |E z (0)/B z (0)| = 0.35, 0.3 and 0.8 respectively. About twenty matching points were sufficient to obtain stable roots. However the field amplitudes were not stable.

We also considered a non-circular base curve of the Pascal's limaçon type (figure 1), parametrized by

r b = a (1 -ε cos s) , ϕ b = s -ε ′ sin s , s ∈ [0, 2π] . ( 21
)
The inward bump is expected to give a stronger E z (0). Indeed when a large part of the base curve is close to the origin, the boundary condition inhibits E z (0). This is avoided with the bump. For ε = ε ′ = 0.3, q = 0.7 and v ph = 0.999, two modes are clearly seen at ωa = 2.4 and 6.7. A state at ωa = 3.5 was also seen for odd l sup , but it appeared only as a smooth dip of ln |det{C nl }| for even l sup . Finally we considered a base curve with two inside bumps and a rotational symmetry of order 2, paramerized by

r b = a (1 -ε cos s) , ϕ b = (s -ε ′ sin s)/2 , s ∈ [0, 4π] . (22)
The modes are either even or odd under this symmetry. Odd modes are not interesting for particle acceleration since E z (0) = 0. Even modes have the property E x (0) = E y (0) = 0, which guaranties that electrons on the zaxis are not deflected and emit no synchrotron radiation. For ε = ε ′ = 0.1, q = 0.7 and v ph = 0.999, even modes are seen at ωa = 4.0 and 6.3, with the different l truncations [l inf , l sup ] = [-4, 8], [-8, 8] and [-10, 10]. With more pronounced bumps (ε = ε ′ = 0.2 or 0.3), no stable root was found.

Limits of the PWE-PM method

We have encountered difficulties when the deformation of the wave guide relative to a straight cylindrical tube becomes too strong. This may be attributed not to the point matching method but to a lack of convergence of the partial wave expansion at r larger than a critical radius r c . Let us indeed consider the expansion [START_REF] Petit | [END_REF], in the scalar field case, and suppose that it converges at some point of polar coordinates (r 0 , ϕ 0 ). Then it converges at all points (r, ϕ) with r < r 0 . It implies that the domain of convergence is bounded by a cylinder.

Demonstration : The convergence of ( 8) at (r 0 , ϕ 0 ) implies that |a l ψ l,k(l) (r 0 )| is bounded by some number M. Using (7), we have

|a l ψ l,k(l) (r)| ≤ M I |l| (κ l r) I |l| (κ l r 0 ) (23)
If r < r 0 , the expansion ( 8) is bound in modulus by a convergent geometric series in |l|, therefore converges at l → ±∞.

If the base curve has a singularity, e.g., an anglar point of polar coordinate (r 1 , ϕ 1 ), the field is non-analytic at this point, therefore r c ≤ r 1 . On the other hand, for a strong enough deformation relative to the straight cylinder, part of the wave guide may be outside the cylinder of convergence, even if the base curve is analytic. A similar situation occurs when one uses the Rayleigh expansion for the wave diffracted by a grating : this expansion does not converge down to the bottom of the grooves if these ones are too deep [START_REF] Petit | [END_REF]9,10].

Conclusion

We have seen that an increase in the strength of the helical bending (increase of ε and q in the case of a de-centered circular basis) leads to a decrease of the dimensionless frequencies ωa of the low modes at imposed phase velocity v ph < c (a characterizes the transverse size of the guide). One can also expect that it leads to a larger accelerating field E z (0), at imposed power of the wave or imposed maximum field on the wall. We have indeed, in order of magnitude

E z /E T ∼ p T /p z , p T ∼ 1/a , p z ∼ ω → E z /E T ∼ 1/(ωa) . ( 25
)
Inward bumps of the base curve also favor large E z (0). On the other hand, the PWE-PM method is limited to weak helical bending and sufficiently "round" base curves. Here we propose two numerical methods to find the propagating modes for larger bending and base curves with more pronounced structures.

The finite element method. This is the standard way to solve numerically the Maxwell equations. In the rotating frame these equations take the 2-dimensional form at fixed ω and P :

iωB X = +∂ Y E z -iP E Y -qE X + q(X∂ Y -Y ∂ X )E Y , iωB Y = -∂ X E z + iP E X -qE Y -q(X∂ Y -Y ∂ X )E X , iωB Z = +∂ X E Y -∂ Y E X , -iP B Z = ∂ X B X + ∂ Y B Y + q(Y ∂ X -X∂ Y )B z , (26)
and analogous equations obtained with the duality transformation B → E, E → -B.

Method of virtual currents. Instead of expanding the field in cylindrical partial waves, one does it with eigenfunctions of P helic (equation 5)which are regular inside the waveguide but singular at some outside places. A candidate is the retarded or advanced field produced by the charge-current distribution

j µ (t, r, z) =   ρ j T j z   = e iP z-iωt δ (2) r -R[qz]d   1 qv ph ẑ ∧ r v ph   , (27)
with v ph = ω/P . R[qz] is the rotation matrix of (1). The x and y components of any vector are gathered in a bi-vector. The point (x, y, z) = (R[qz]d, z) draws an helix of path Λ = 2π/q intersecting the z = 0 plane at point (d, 0). ρ(t, x, y, z) is a periodic distribution of charges running along this helix with the longitudinal velocity v z = v ph and j is the associated current. j µ is periodic in time with frequency ω and is an eigenfunction of P helic with eigenvalue P , therefore the retarded or advanced field has the same property. Let us denote, for instance, by F (d) = {E (d) , B (d) } the retarded field. Different d's give different independent solutions. To take into account the fact that there are two polarisations, we also consider the dual field F(d) = {B (d) , -E (d) }, which would be produced by a distribution of magnetic charge and current. In the z = 0 plane these fields have a singularity in 1/|r -d| at r = d ; they are Green's functions of (26). Choosing N singular points {d 1 , d 2 , • • • d N , } external to the base curve, we may approximate the field of a mode by

F = N n=1 a n F (dn) + b n F(dn) , (28)
The coefficients a n and b n should be choosen to satisfy the boundary conditions on N matching points r 1 , ...r N of the base curve. Near a sharp inward bump of the base curve the field is strong and, by analytic continuation, one probably meets a singularity not far behind the curve, analogous to a point-like electrical image in electrostatics. This singularity may be simulated by placing one or a few d n 's on it or near to it. It amounts to replace the continuous surface charge and currents on the bump by charges and currents on virtual helical wires.

Whichever method will be used, it will be usefull to compare it with the PWE-MP method in the case of weak helical deformation, as a check.

  2.2). For a base without discrete rotational symmetry, L = [l inf , l sup ] with l sup = l inf + N -1. The choice of r n 's and L should be guided by the following conditions :

Figure 2 :

 2 Figure 2: 3-dimensional representation of the scalar field Ψ(x, y) = Φ(0, x, y, 0) of the ωa = 5.0 mode, for a circular base of radius a, off-centred by ε = 0.4 from the axis, and q = 1. Left: real part; right: imaginary part.

  8 and 5.0. They remain the same for various L sets like [-6, +8], [-5, +9] and [-3, +11]. The TM/TE character, measured by |E z (0)/B z (0)|, are 0.3, 0.25 and 0.4. and 5.4 respectively. The ωa = 5.0 mode seems therefore more appropriate for electron acceleration. It has a figure of merit f.o.m. = |E z (0)|/[max(|E|) on the wall] ≃ 1/8 (20) and a group velocity v g = dω/dP = 0.68. Plots of the fields of this mode can be seen in [7]. Let us mention that |B r | and |B ϕ | look very similar to |E ϕ | and |E r | respectively, which is reminiscent of a plane wave.

  for any r. At large |l|, we have κ l ≃ |l|q -s l P with s l = sign(l) and the asymptotic form of the modified Bessel functions gives I |l| (κ l r)/I |l| (κ l r 0 ) ≃ exp[(|l|q -s l P ) (r -r 0 )]. Then |a l ψ l,k(l) (r)| ≤ M exp[s l P (r 0 -r)] e |l|q(r-r 0 ) at large l .

Presented at the VIII th International Symposium Radiation from Relativistic Electrons in Periodic Structures (RREPS-09), Zvenigorod, Russia, Sept. 7-11,


X Artru
email: x.artru@ipnl.in2p3.fr

RECURSIVE FRAGMENTATION MODEL WITH QUARK SPIN. APPLICATION TO QUARK POLARIMETRY 1

An elementary recursive model accounting for the quark spin in the fragmentation of a quark into mesons is presented. The quark spin degree of freedom is represented by a two-components spinor. Spin one meson can be included. The model produces Collins effect and jet handedness. The influence of the initial quark polarisation decays exponentially with the rank of the meson, at different rates for longitudinal and transverse polarisations.

Introduction

Present Monte-Carlo event generators of quark and gluon jets do not include the parton spin degree of freedom, therefore do not generate the Collins [1] and jet handedness [2] effects. These are azimuthal asymmetries bearing on one, two or three hadrons, which can serve as quark polarimeters. However the asymmetries may strongly depend, in magnitude and sign, on the quark and hadron flavors and on the transverse momenta p T and scaled longitudinal momenta z of these hadrons. Therefore a good knowledge of this dependence is needed for parton polarimetry. Due to the large number of kinematical variables, a hadronisation model which takes spin into account is urgently needed as a guide.

The semi-classical Lund 3 P 0 mechanism [3], grafted on the string model, can generate a Collins effect [4], but not jet-handedness. Here we propose a fully quantum model of spinning quark fragmentation, based on the multiperipheral model. It reproduces the results of the 3 P 0 mechanism and also contains the jet-handedness effect.

2 Some recalls about quark fragmentation Figure 1 describes the creation of a quark "q 0 " and an antiquark "q -1 " in e + e -annihilation or W ± decay, followed by the hadronisation,

q 0 + q-1 → h 1 + h 2 ... + h N . (1)
Looking from rigth to left, one sees it as the recursive process (see [5] and ref. 4 of [6]), q N ≡ q -1 is a "quark propagating backward in time" and k N ≡ -k(q -1 ).

q 0 ≡ q 0 → h 1 + q 1 q 1 → h 2 + q 2 • • • q N -1 → h N + q N 4-momenta : k 0 = p 1 + k 1 , k 1 = p 2 + k 2 , • • • k N -1 = p N + k N .

Kinematical notations :

k 0 = k(q 0 ) and k(q -1 ) are in the +ẑ and -ẑ directions respectively. For a quark, t n ≡ k nT . For a 4-vector, a ± = a 0 ± a z and a T = (a x , a y ). We denote by a tilde the dual transverse vector ãT ≡ ẑ × a T = (-a y , a x ).

In Monte-Carlo simulations, the k n are generated according to the splitting distribution

dW ( q n-1 → h n + q n ) = f n (ζ n , t 2 n-1 , t 2 n , p 2 nT , ) dζ n d 2 t n , ζ n ≡ p + n /k + n-1 .
In particular the symmetric Lund splitting function [3],

f n ∝ ζ a n-1 -an-1 n (1 -ζ an ) exp -b (m 2 n + p 2 nT )/ζ n , (3)
inspired by the string model, fulfills the requirement of forward-backward equivalence.

On can also consider [6] the upper part of Fig. 1 as a multiperipheral [7] diagram with the Feynman amplitude

M q 0 +q -1 →h 1 ...+h N = v( k -1 , S -1 ) Γ q N ,h N ,q N-1 (k N , k N -1 ) ∆ q N-1 (k N -1 ) • • • • • • ∆ q 2 (k 2 ) Γ q 2 ,h 2 ,q 1 (k 2 , k 1 ) ∆ q 1 (k 1 ) Γ q 1 ,h 1 ,q 0 (k 1 , k 0 ) u(k 0 , S 0 ) . (4)
S 0 and S -1 are the polarisation vectors of the intial quark and antiquark. S 2 = 1, S z = helicity, S T = transversity. Γ and ∆ are vertex functions and propagators which depend on the quark momenta and flavors. Note that Fig. 1 is a loop diagram : k 0 is an integration variable, therefore the "jet axis" is not really defined. Furthermore, in Z 0 or γ * decay, the spins q 0 and q-1 are entangled so that one cannot define S 0 and S -1 separately.

Collins and jet-handedness effects. Let us first assume that the jet axis (quark direction) is well determined :

-the Collins effect [1], in q → h+X, is an asymmetry in sin[ϕ(S)-ϕ(h)] for a transversely polarized quark. The fragmentation function reads

F (z, p T ; S T ) = F 0 (z, p 2 T ) (1 + A T S T .p T /|p T |) (p T ≡ ẑ × p T ) . (5)
A T = A T (z, p 2 T ) ∈ [-1, +1
] is the Collins analysing power. -jet handedness [2], in q → h + h ′ + X, is an asymmetry in sin[ϕ(h) -ϕ(h ′ )] proportional to the quark helicity. The 2-particle longitudinaly polarised fragmentation function is

F (z, p T , z ′ , p ′ T ; S z ) = F 0 (z, p 2 T , z ′ , p ′ T 2 , p T • p ′ T ) 1 + A L S z pT .p ′ T |p T • p ′ T | . ( 6
)
A L = A L (z, p 2 T , z ′ , p ′ T 2 , p T .p ′ T ) ∈ [-1, +1
] is the handedness analysing power. p1T • p ′ T is the same as ẑ • (p T × p ′ T ). If the jet axis is not well determined, an additional fast hadron, h ′ or h ′′ is needed. The z axis is taken along P = (p+p ′ ) (Collins) or P = (p+p ′ +p ′′ ) (handedness). In this way, we define the 2-particle relative Collins effect (also called interference fragmentation) and the three-particle jet handedness, which corresponds to the original definition of [2].

The Lund 3 P 0 mechanism [3]. Figure 2 depicts the decay of the initial massive string accompagnied with the creation of a q q pairs. Forgetting transverse oscillations of the initial string, the transverse hadron momenta come from the internal orbital motions of the pairs. After a tunnel effect the q and q of a pair become on-shell and their relative position r ≡ r(q)r(q) is along -ẑ. The pair is assumed to be in the 3 P 0 state, which has the vacuum quantum number. The relative momentum k ≡ k(q) = -k(q) and the orbital angular momentum L = r × k are such that ẑ • [k T × L] < 0. In the 3 P 0 state s q = s q = -L/2 . As a result, the transverse spins of q and q are correlated to their transverse momenta :

ẑ • [k T (q) × s q ] > 0 , ẑ • [k T (q) × s q] < 0 . (7)
The correlation can be transmitted to a baryon. Then ẑ • [p T × s B ] has the sign of s q .s B . This can explain transverse spin asymmetries in hyperon production [3].

Application to the Collins effect [4]. In Fig. 2, q 0 is polarised along the direction +ŷ toward the reader and h 1 is a pseudoscalar meson, for which s(q 0 ) =s(q 1 ) . Then q 1 and q1 are polarised along -ŷ and, according to (7), k T (q 1 ) = p T (h 1 ) is in the +x direction. This provides a model for the Collins effect. Fig. 2 also indicates that, for a sequence of pseudoscalar mesons, the Collins asymmetries are of alternate sides. Besides, q n-1 and qn go on the same side, which enhances the asymmetry. It may explain why π - from u-quarks have a strong Collins analysing power. Note that this effect also enhances p 2 T , independently of the q 0 polarisation.

A simplified multiperipheral quark model

In Eq.( 4), let us replace Dirac spinors by Pauli spinors. A minimal model, restricted to the direct emission of pseudoscalar mesons, is built with the following prescriptions : 1) replace u(k 0 , S 0 ) and v(k -1 , S -1 ) ≡ -ū(k q-1 , -S q-1 ) γ 5 by the Pauli spinors χ(S 0 ) and

-χ † (-S -1 ) σ z , 2) assume no momentum dependence of Γ , 3) replace γ 5 by σ z , 4) replace the usual pole (k 2 -m 2 q ) -1 of ∆ q (k) by the (k L , k T ) separable form D q (k) = g q (k + k -) exp(-Bt 2 /2) , (8) 5
) replace the usual numerator m q + γ.k by µ q (k + k -, t 2 ) + iσ. t . These prescriptions respect the invariance under the following transformations :

-rotation about the z-axis, -Lorentz transformations along the z-axis (longitudinal boost), -mirror reflection about any plane containing the z-axis (parity), -forward-backward equivalence.

The jet axis being fixed, full Lorentz invariance is not required, whence the separate dependences of D q and µ q in k + k -and t 2 . In item 5), µ+iσ. t is reminiscent of the mesonbaryon scattering amplitude f (s, t)+ig(s, t) σ.(p×p ′ ). Single-spin effects are obtained for ℑm µ = 0. The choice of putting the spin dependence in the propagators rather than in the vertices is inspired by the 3 P 0 mechanism : in both models the polarisation germinates in the quark line between two hadrons. For a fast investigation of the model, we make the further approximations : -Neglect the influence of the antiquark flavor and polarisation in the quark fragmentation region. This is allowed at large invariant q 0 + q-1 mass.

-Discard the interference diagrams. For a given final state, the rank ordering of hadrons in the multiperipheral diagram is not unique and differently ordered diagrams can interfere. This interference (and the resulting Bose-Einstein correlations) will be neglected.

-Disentangle k ± and k T . We will assume that µ q (k + k -, t 2 ) is constant or a function of t 2 only. Thus we have no more "dynamical" correlation between longitudinal and transverse momenta. However there remains a "kinematical" correlation coming from the mass shell constraint

(k n-1 -k n ) 2 ≡ (k + n-1 -k + n )(k - n-1 -k - n ) -(t n-1 -t n ) 2 = m 2 (h n ) . (9)
In the following we will ignore the (t n-1t n ) 2 term. This approximation is drastic for pion emission because t 2 > m 2 π . We only use it here for a qualitative investigation of the spin effects allowed by the multiperipheral model. Thanks to it, the t's become fully decoupled from the k ± n and kinematically decorrelated between themselves. They remain correlated only via the quark spin.

p T -distibutions in the quark fragmentation region. With the above approximations we can treat the process (2), at least in p T -space, like a cascade decay of unstable particles, which has no constraint coming from the future. The joint p T distibution of the n first mesons is proportional to

I(p 1T , p 2T , ...p nT ) = exp(-B t 2 1 -B t 2 2 ... -B t 2 n ) Tr M 12...n 1 + S 0 .σ 2 M † 12...n , (10) with
M 12...n = M n • • • M 2 M 1 , M r = (µ r + iσ. tr ) σ z . (11)

Applications to azimuthal asymmetries

In this section we will calculate azimuthal asymmetries for particles of definite ranks. For comparison with experiments, one should mix the contributions of different rank assignments. For simplicity we take a unique and constant µ for all quark flavors.

First-rank Collins effect. Applying (10)(11) for n = 1 gives

I(p 1T ) = exp(-Bt 2 1 ) |µ| 2 + t 2 1 -2 ℑm(µ) t1 .S , (12)
with t 1 = -p 1T . For complex µ one has a Collins asymmetry (cf Eq.5) with

A T = 2 ℑm(µ) |p 1,T | |µ| 2 + p 2 1,T ∈ [-1, +1] . (13)
If ℑm(µ) > 0 it has the same sign as predicted by the 3 P 0 mechanism.

Joint p T spectrum of h 1 and h 2 . Applying (10)(11) for n = 2 one obtains

I(p 1T , p 2T ) = exp(-Bt 2 1 -Bt 2 2 ) { (|µ| 2 + t 2 1 ) (|µ| 2 + t 2 2 ) -4t 1 .t 2 ℑm 2 (µ) + 2 ℑm(µ) S. t1 (2 t 1 .t 2 -|µ| 2 -t 2 2 ) + 2 ℑm(µ) S. t2 (|µ| 2 -t 2 1 ) -2 ℑm(µ 2 ) S.(t 1 × t 2 ) } , (14)
with

t 1 = -p 1T , t 2 = -(p 1T + p 2T ) and S • (t 1 × t 2 ) = S z p1T • p 2T .
The last line contains jet handedness (cf Eq.6), of analysing power

A L = -2 ℑm(µ 2 ) |p 1T × p 2T | (|µ| 2 + t 2 1 ) (|µ| 2 + t 2 2 ) -4t 1 .t 2 ℑm 2 (µ) ∈ [-1, +1] . ( 15
)
The second line contains the Collins asymmetry of h 1 . Both 2 nd and 3 rd lines contribute to the h 2 one after integration over t 1 , and to the relative 2-particle Collins asymmetry, which bears on

r 12 = z 2 p 1T -z 1 p 2T z 1 + z 2 = z 1 z 1 + z 2 t 2 -t 1 . ( 16
)
Note that Collins and jet-handedness asymmetries are not maximum for the same value of arg(µ). This is related to the positivity [8] constraint

A 2 L (p 1T , p 2T ) + A 2 T (p 1T , p 2T ) ≤ 1 . ( 17
)

Evolution of the polarisation of the cascading quark

Let us first assume that t 1 , t 2 , ... t n are fixed and consider the spin density matrix ρ n = (1 + S n .σ)/2 of q n at the (n + 1) th vertex :

ρ n = R n / Tr{R n } , R n = M 12...n 1 + S 0 .σ 2 M † 12...n . ( 18
)
If ρ 0 is a pure state (det ρ 0 = 0), then ρ n is also a pure state ; no information is lost.

Let us now integrate over t 1 , t 2 , ... t n (equivalently over p 1T , ...p nT ). It leads to a loss of information. The spin density matrix of q n becomes ρn = Rn / Tr{ Rn } ,

Rn = d 2 t 1 ... d 2 t n M 12...n 1 + S 0 .σ 2 M † 12...n . (19)
R n and Rn obey the recursion relations

R n = M n R n-1 M † n , Rn = d 2 t n M n Rn-1 M † n . ( 20
)
At t's, the left equation gives (setting µ = µ ′ + iµ ′′ ) :

S n = 1 C 2µ ′′ tn + R[ẑ, ϕ n ]   t 2 -|µ| 2 0 -2|t| µ ′ 0 -|µ| 2 -t 2 0 2|t| µ ′ 0 |µ| 2 -t 2   R[ẑ, -ϕ n ] S n-1 , (21)
with (12), where we replace {S, t 1 } by {S n-1 , t n }, and (21), we generate the successive transverse momenta with the Monte-Carlo method. From (21) we learn :

C = Tr{R n } = |µ| 2 + t 2 n -2µ ′′ tn • S n-1 . The rotation R[ẑ, ϕ n ] about ẑ brings x along t n . Iterating
-if ℑm µ = 0, the inhomogeneous term in µ ′′ tn is a source (or sink) of transverse polarisation : one can have S nT = 0 even with S n-1 = 0.

-helicity is partly converted into transversity along t n and vice-versa.

The last fact explains the mechanism of jet handedness in this model : first, the helicity S z0 is partly converted into S 1T parallel to p 1T , then S 1T produces a Collins asymmetry for h 2 in the plane perpendicular to p 1T .

Let us now consider the t-integrated density matrix. The right equation in (20) gives

S n,z = D LL S n-1,z , S n,T = D T T S n-1,T ; D LL , D T T ∈ [-1, +1] , (22)
D LL D T T = d 2 t exp(-Bt 2 ) |µ| 2 -t 2 -|µ| 2 d 2 t exp(-Bt 2 ) (|µ| 2 + t 2 ) . (23)
Indeed, 2D T T = -1 -D LL . This is due to the zero spin of h n (compare with text after Eq.(4.87) of [8]). The negative value of D T T leads to Collins asymmetries of alternate signs, in accordance with the 3 P 0 mechanism. It comes from the σ z vertex for pseudoscalar mesons. For scalar mesons we replace σ z by 1. In this case D T T is positive, q n-1 and qn tend towards opposite sides and the Collins effect is small, except for h 1 . This is also the prediction of the 3 P 0 mechanism.

Inclusion of spin-1 mesons

For a J P C = 1 --vector meson and the associated self-conjugate multiplet, the "minimal" emission vertex written with Pauli matrices is

Γ = G L V * z 1 + G T σ.V * T σ z , ( 25
)
where V is the vector amplitude of the meson normalised to V .V * = 1. It is obtained from the relativistic 4-vector V µ first by a longitudinal boost which brings the hadron at p z = 0, then a transverse boost which brings the hadron at rest. For a J P C = 1 ++ axial meson of amplitude A, the "minimal" emission vertex is

Γ = GT σ.A * T . (26)
It differs by a σ z matrix from the second term of (25). A term of the form GL A * z σ z with constant GL is not allowed by the forward-backward equivalence.

Let us treat the case where the 1 st -rank particle is a ρ + meson and fix the momenta p(π + ) and p(π 0 ) of the decay pions. Then V µ ∝ p(π + )-p(π 0 ), which is real, corresponding to a linear polarisation. Replacing the σ z coupling of ( 11) by ( 25) we obtain

I(p T , V ) = exp(-Bt 2 ) |G T | 2 × { (|α| 2 V 2 z + V 2 T ) (|µ| 2 + t 2 ) -4V T .t V z ℑm(α) ℑm(µ) + 2 ℑm(µ) |α| 2 V 2 z S. t + 2 ℑm(α) (|µ| 2 + t 2 ) V z V T . S + 2 ℑm(µ) (V T . t V T .S + V T .t V T . S) + 4 ℜe(α) ℑm(µ) S z V z V T . t } , (27)
with t ≡ t 1 = -p T (ρ + ) and α ≡ G L /G T . Let us comment this formula :

• The 2 nd line is for unpolarised quark. It gives some tensor polarisation.

• The 3 rd line is a Collins effect for the ρ + as a whole, opposite to the pion one (compare with (12) and only for longitudinal linear polarisation, in accordance with the 3 P 0 mechanism. For V 2 z = 1/3 (unpolarized ρ + ) and α = 1 one recovers the Czyzewski prediction [9] A T (leading ρ)/A T (leading π) = -1/3. • The 4 th line gives an oblique polarisation in the plane perpendicular to S T corresponding to ĥ1 or h 1LT of [10,11]. After ρ + decay, it becomes a relative π + -π 0 Collins effect.

• The 5 th line is a new type of asymmetry, in sin[2ϕ(V ) -ϕ(t) -ϕ(S)].

• The last line also gives an oblique polarisation, but in the plane perpendicular to p T (ρ + ). After ρ + decay, it becomes jet-handedness. Indeed, ignoring an effect of transverse boost, we have

V T ∝ p T (π + ) -p T (π 0 ), therefore V T . t ∝ -p T (π + ) × p T (π 0 ).

Conclusion

For the direct fragmentation of a transversely polarised quark into pseudo-scalar mesons, the model we have presented has essentially one free complex parameter µ and reproduces the results of the semi-classical 3 P 0 mechanism : large asymmetry for the 2 nd -rank meson, Collins asymmetries of alternate sides for the subsequent mesons. In addition, it possesses a jet-handedness asymmetry, generated in two steps : partial transformation of helicity into transversity, then Collins effect.

We have also considered the inclusion of spin-1 mesons. When longitudinally polarised, a leading vector meson has a Collins asymmetry opposite to that of a pseudoscalar, as also expected from the 3 P 0 mechanism. The decay pions of a ρ meson exhibit a relative Collins effects as well as jet-handedness. These effects are associated to oblique linear polarisations of the ρ meson. The fact that two pions coming from a ρ show the same spin effects as two successive "direct" pions is reminiscent of duality.

Even for unpolarised initial quarks, the spin degree of freedom of the cascading quark has to be considered. It enhances the p 2 T of the pseudoscalar mesons compared to scalar and longitudinal vector mesons.

A next task for building a realistic Monte-Carlo generator with quark spin is to take into account the (t n-1t n ) 2 term in (9). One must also be aware that there exist other mechanisms of spin asymmetries in jets. For example the Collins effect can be generated by the interference between direct emission and the emission via a resonance [12].

( 2 )Figure 1 .

 21 Figure 1. Electroweak boson → q q → mesons.

Figure 2 .

 2 Figure 2. String decaying into pseudoscalar mesons.

  Analytical values : D LL = (ξ -1)/(ξ + 1) and D T T = -ξ/(ξ + 1) with ξ = B|µ| 2 . The geometrical decays of |S n,z | and |S n,T | along the quark chain occur at different speeds. They are similar to the decays of charge and strangeness correlations. D LL and D T T saturate a Soffer-type [8] positivity condition 2|D T T | ≤ 1 + D LL .


Manuel Caroli

Pedro M M De Castro

Sébastien Loriot
email: sebastien.loriot@sophia.inria.fr

Olivier Rouiller

Monique Teillaud
email: monique.teillaud@sophia.inria.fr

Camille Wormser
email: camille.wormser@inf.ethz.ch

{manuel Caroli

Pedro Machado

Robust and Ecient Delaunay triangulations of points on or close to a sphere

Keywords: Computational Geometry, Convex hull, Delaunay Triangulation, Voronoi Diagram, Sphere, Space of Circles, Exact Geometric Computing, cgal Géométrie algorithmique, Enveloppe convexe, Triangulation

We propose two ways to compute the Delaunay triangulation of points on a sphere, or of rounded points close to a sphere, both based on the classic incremental algorithm initially designed for the plane. We use the so-called space of circles as mathematical background for this work. We present a fully robust implementation built upon existing generic algorithms provided by the cgal library. The eciency of the implementation is established by benchmarks.

Nous proposons deux façons de calculer la triangulation de Delaunay d'un ensemble de points qui appartiennent soit à la sphère, soit à son voisinage. Ces deux méthodes reposent sur l'algorithme incrémental classique, tel qu'il a été créé à l'origine pour calculer les triangulations de Delaunay planaires. Le cadre mathématique classique justiant cette approche est rappelé, à l'aide de l'espace des cercles. Ces deux approches ont été implantées de façon robuste en s'appuyant sur les algorithmes génériques fournis par la bibliothèque CGAL. Des tests comparatifs montrent l'ecacité de nos implantations sur des jeux de données de taille variée.

Introduction

The cgal project [cga] provides users with a public discussion mailing list, where they are invited to post questions and express their needs. There are recurring requests for a package computing the Delaunay triangulation or its dual, the Voronoi diagram, of points on a sphere. This is useful in many domains such as geology, geographic information systems or structural molecular biology, to name a few. An easy and standard solution to the problem of computing such a Delaunay triangulation consists in constructing the 3D convex hull of the points: They are equivalent [START_REF] Brown | Geometric transforms for fast geometric algorithms[END_REF][START_REF] Sugihara | Laguerre Voronoi diagram on the sphere[END_REF]. The convex hull is one of the most popular structures in computational geometry [START_REF] Mark De Berg | Computational Geometry: Algorithms and Applications[END_REF][START_REF] Boissonnat | Algorithmic Geometry[END_REF]; optimal algorithms and ecient implementations are available [hul, qhu].

Another fruitful way to compute Delaunay on a sphere consists of reworking known algorithms designed for computing triangulations in R 2 . Renka adapts the distance in the plane to a geodesic distance on a sphere and triangulates points on a sphere [START_REF] Renka | Algorithm 772: STRIPACK: Delaunay triangulation and Voronoi diagram on the surface of a sphere[END_REF] through the well-known ipping algorithm for Delaunay triangulations in R 2 [START_REF] Lawson | Software for C 1 surface interpolation[END_REF]. As a by-product of their algorithm for arrangements of circular arcs, Fogel et al. can compute Voronoi diagrams of points lying exactly on the sphere [FSH08, FS, BFH + 09b, BFH + 09a]. Using two inversions allows Na et al. to reduce the computation of a Voronoi diagram of sites on a sphere to computing two Voronoi diagrams in R 2 [NLC02], but no implementation is available. Note that this method assumes that data points are lying exactly on a sphere.

As we are motivated by applications, we take practical issues into account carefully. While data points lying exactly on the sphere can be provided either by using Cartesian coordinates represented by a number type capable of handling algebraic numbers exactly, or by using spherical coordinates, in practice data-sets in Cartesian coordinates with double precision are most common. In this setting, the data consists of rounded points that do not exactly lie on the sphere, but close to it.

In Section 4, we propose two dierent ways to handle such rounded data. Both approaches adapt the well-known incremental algorithm [START_REF] Bowyer | Computing Dirichlet tessellations[END_REF] to the case of points on, or close to the sphere. It is important to notice that, even though data points are rounded, we follow the exact geometric computation paradigm pioneered by C. K. Yap [START_REF] Yap | The exact computation paradigm[END_REF]. Indeed, it is now well understood that simply relying on oating point arithmetic for algorithms of this type is bound to fail (see [KMP + 08] for instance).

The rst approach (Section 4.1) considers as input the projections of the rounded-data points onto the sphere. Their coordinates are algebraic numbers of degree two. The approach computes the Delaunay triangulation of these points exactly lying on the sphere.

The second approach (Section 4.2) considers circles on the sphere as input. The radius of a circle (which can alternatively be seen as a weighted point) depends on the distance of the corresponding point to the sphere. The approach computes the weighted Delaunay triangulation of these circles on the sphere, also known as the regular triangulation, which is the dual of the Laguerre Voronoi diagram on the sphere [START_REF] Sugihara | Laguerre Voronoi diagram on the sphere[END_REF] and the convex hull of the rounded-data points. These interpretations of rounded data are supported by the space of circles [Ber87, DMT92] (Section 3).

The same algorithm being used for both approaches, only the interpretation of the data and the way in which basic computations on them are performed modify the behavior and possibly the result of the algorithm: Some circles may be hidden in a regular triangulation. We show in Section 4.2 that this can happen only if the rounded points are extremely close to one another. Thus, under some sampling conditions, the two approaches compute the same triangulation (except maybe in degenerate cases for which the considered triangulations are known not to be uniquely dened), and in general, all points appear in the result of both approaches.

In Section 5, we show that the genericity of cgal [START_REF] Fogel | Generic programming and the CGAL library[END_REF] enables the reuse of a large part of the existing package for triangulations in R 2 . We present experimental results on very large data-sets, showing the eciency of our approaches. We compare our code to software designed for computing Delaunay triangulations on the sphere, and to convex-hull software [HS09, PT09b, hul, sug, qhu, Ren97, FS]. The performance, robustness, and scalability of our approaches express their added value.

Denitions and Notation

Let us rst recall the denition of the regular triangulation in R 2 , also known as weighted Delaunay triangulation. A circle c with center p ∈ R 2 and squared radius r 2 is considered equivalently as a weighted point and is denoted by c = (p, r 2 ). The power product of c = (p, r 2 ) and c = (p , r 2 ) is dened as pow(c, c ) = pp 2 -r 2 -r 2 , where pp denotes the Euclidean distance between p and p . Circles c and c are orthogonal i pow(c, c ) = 0. If pow(c, c ) > 0 (i.e., the disks dened by c and c do not intersect, or the circles intersect with an angle strictly smaller than π

2 ), we say that c and c are suborthogonal. If pow(c, c ) < 0, then we say that c and c are superorthogonal (see Figure 1). Three circles whose centers are not collinear have a unique common orthogonal circle. Let S be a set of circles. Given three circles of S, c i = (p i , r 2 i ), i ∈ I ⊂ N and |I| = 3, whose centers are not collinear, let T I be the triangle whose vertices are the three centers p i , i ∈ I. We dene the orthogonal circle of T I as the circle that is orthogonal to the three circles c i , i ∈ I. T I is said to be regular if for any circle c ∈ S, the orthogonal circle of T I and c are not superorthogonal. A regular triangulation RT (S) is a partition of the

INRIA π/2 < π/2 > π/2
Figure 1: From left to right: orthogonal (pow(s 0 , s 1 ) = 0), suborthogonal (pow(s 0 , s 1 ) > 0), and superorthogonal (pow(s 0 , s 1 ) < 0) circles in R 2 . If all radii are equal, then the power test reduces to testing whether a point lies inside, outside, or on the circle passing through three points; the regular triangulation of the circles is the Delaunay triangulation DT of their centers.

More background can be found in [START_REF] Aurenhammer | Power diagrams: properties, algorithms and applications[END_REF]. We refer the reader to standard textbooks for algorithms computing Delaunay and regular triangulations [START_REF] Mark De Berg | Computational Geometry: Algorithms and Applications[END_REF][START_REF] Boissonnat | Algorithmic Geometry[END_REF].

This denition generalizes in a natural manner to the case of circles lying on a sphere S in R 3 : Angles between circles are measured on the sphere, triangles are drawn on the sphere, their edges being arcs of great circles. As can be seen in the next section, the space of circles provides a geometric presentation showing without any computation that the regular triangulation on S is a convex hull in R 3 [START_REF] Sugihara | Laguerre Voronoi diagram on the sphere[END_REF].

In the sequel, we assume that S is given by its center, having rational coordinates (we take the origin O without loss of generality), and a rational squared radius R 2 . This is also how spheres are represented in cgal. 1 1 We mention rational numbers to simplify the presentation. cgal allows more general number types that provide eld operations: +, -, ×, /.

RR n°7004

3 Space of Circles Computational geometers are familiar with the classic idea of lifting up sites from the Euclidean plane onto the unit paraboloid Π in R 3 [START_REF] Aurenhammer | Voronoi diagrams: A survey of a fundamental geometric data structure[END_REF]. We quickly recall the notion of space of circles here and refer to the literature for a more detailed presentation [START_REF] Devillers | The space of spheres, a geometric tool to unify duality results on Voronoi diagrams[END_REF]. In this lifting, points of R 3 are viewed as circles of R 2 in the space of circles: A circle c = (p, r 2 ) in R 2 is mapped by π to the point π(c) = (x p , y p , x 2 p + y 2 p -r 2 ) ∈ R 3 . A point of R 3 lying respectively outside, inside, or on the paraboloid Π represents a circle with respectively positive, imaginary, or null radius. The circle c in R 2 corresponding to a point π(c) of R 3 outside Π is obtained as the projection onto R 2 of the intersection between Π and the cone formed by lines through π(c) that are tangent to Π; this intersection is also the intersection of the polar plane P (c) of π(c) with respect to the quadric Π.

Points lying respectively on, above, below P (c) correspond to circles in R 2 that are respectively orthogonal, suborthogonal, superorthogonal to c.

The predicate pow(c, c ) introduced above is thus equivalent to the orientation predicate in R 3 that tests whether the point π(c ) lies on, above or below the plane P (c). If c is the common orthogonal circle to three input circles c 1 , c 2 , and c 3 (where

c i = (p i , r 2 i ) for each i), then pow(c, c ) is the orientation predicate of the four points π(c 1 ), π(c 2 ), π(c 3 ), π(c ) of R 3 . It can be expressed as sign ˛1 1 1 1 xp 1 xp 2 xp 3 x p yp 1 yp 2 yp 3 y p zp 1 zp 2 zp 3 z p ˛, (1)
where z p i = x 2 p i + y 2 p i -r 2 i for each i and z 2 p = x 2 p + y 2 p -r 2 . It allows to relate Delaunay or regular triangulations in R 2 and convex hulls in R 3 [Aur91], while Voronoi diagrams in R 2 are related to upper envelopes of planes in R 3 .

Up to a projective transformation, a sphere in R 3 can be used for the lifting instead of the usual paraboloid [START_REF] Berger | The space of spheres[END_REF]. In this representation, common in geometry, the sphere has a pole 2 and can be identied to the Euclidean plane R 2 . What we are interested in this paper is the space of circles drawn on the sphere S itself, without any pole. This space of circles has a nice relation to the de Sitter space in Minkowskian geometry [START_REF] Coxeter | A geometrical background for de Sitter's world[END_REF]. We can still construct the circle c on S that is associated to a point p = π S (c) of R 3 as the intersection between S and the polar plane P S (p) of p with respect to the quadric S (Figure 3). Its center is the projection of p onto S and as above, imaginary radii are possible. 3 So, in the determinant in (1), x p i , y p i , and z p i (respectively x p , y p , z p ) are precisely the coordinates of the points p i = π S (c i ) (respectively p = π S (p)). This will be extensively used in Sec-2 See the nice treatment of innity in [START_REF] Berger | The space of spheres[END_REF].

p 2 = π S (c 2 ) c 2 Figure 3: c 1 is suborthogonal to c, c 2 is superorthogonal to c.
tion 4. Again, we remark that Delaunay and regular triangulations on S relate to convex hulls in 3D. Interestingly, rather than using a convex hull algorithm to obtain the Delaunay or regular triangulation on the surface as usually done for R 2 [Aur91], we will do the converse in the next section.

Algorithm

The incremental algorithm for computing a regular triangulation of circles on the sphere S is a direct adaptation of the algorithm in R 2 [START_REF] Bowyer | Computing Dirichlet tessellations[END_REF]. Assume that RT i-1 = RT ({c j ∈ S, j = 1, . . . , i -1}) has been computed. 4 The insertion of c i = (p i , r 2 i ) works as follows: • locate p i (i.e., nd the triangle t containing p i ),

• if t is hiding p i (i.e., if c i and the orthogonal circle of t are suborthogonal) then stop; p i is not a vertex of RT i . Note that this case never occurs for Delaunay triangulations.

• else (i) nd all triangles whose orthogonal circles are superorthogonal to c i and remove them; this forms a polygonal region that is star-shaped with respect to p i ; 5 (ii) triangulate the polygonal region just created by constructing the triangles formed by the boundary edges of the region and the point p i .

Two main predicates are used by this algorithm:

The orientation predicate allows to check the orientation of three points 4 For the sake of simplicity, we assume that the center O of S lies in the convex hull of the data-set. This is likely to be the case in practical applications. So, we just initialize the triangulation with four dummy points that contain O in their convex hull and can optionally be removed in the end.

5 As previously noted for the edges of triangles, all usual terms referring to segments are transposed to arcs of great circles on the sphere.

p, q, and r on the sphere. (This predicate is used in particular to locate new points.) It is equivalent to computing the side of the plane dened by O, p, and q on which r is lying, i.e., the orientation of O, p, q, and r in R 3 .

The power test introduced in Section 2 boils down to an orientation predicate in R 3 , as seen in Section 3. (This predicate is used to identify the triangles whose orthogonal circles are superorthogonal to each new circle.)

The two approaches briey presented in the introduction fall into the general framework of computing the regular triangulation of circles on the sphere. The next two sections precisely show how these predicates are evaluated in each approach.

4.1

First approach: using points on the sphere

In this approach, input points for the computation are chosen to be the projections on S of the rounded points of the data-set with rational coordinates. The three coordinates of an input point are thus algebraic numbers of degree two lying in the same extension eld of the rationals. In this approach weights, or equivalently radii if circles, are null. The power test consists in this case in answering whether a point s lies inside, outside, 6 or on the circle passing through p, q, and r on the sphere. Following Section 3, this is given by the orientation of p, q, r, and s, since points on the sphere are mapped to themselves by π S .

The diculty comes from the fact that input points have algebraic coordinates. The coordinates of two dierent input points on the sphere are in general lying in dierent extensions. Then the 3D orientation predicate of p, q, r, and s given by (1) is the sign of an expression lying in an algebraic extension of degree 16 over the rationals, of the form

a 1 √ α 1 + a 2 √ α 2 + a 3 √ α 3 + a 4 √ α 4
where all a's and α's are rational. Evaluating this sign in an exact way allows to follow the exact computation framework ensuring the robustness of the algorithm.

Though software packages oer exact operations on general algebraic numbers [cor, led], they are much slower than computing with rational numbers. The sign of the above simple expression can be computed as follows: 1 evaluate the signs of

A 1 = a 1 √ α 1 +a 2 √ α 2 and A 2 = a 3 √ α 3 +a 4 √ α 4 , by comparing a i √ α i with a i+1
√ α i+1 for i = 1, 3, which reduces after squaring to comparing two rational numbers, 2 the result follows if A 1 and A 2 have dierent signs, 3 otherwise, compare A 2 1 with A 2 2 , which is an easier instance of 1. To summarize, the predicate is given by the sign of polynomial expressions in the rational coordinates of the rounded-data points, which can be computed exactly using rational numbers only.

4.2

Second approach: using weighted points

In this approach, the regular triangulation of the weighted points is computed as described above. As in the previous approach, both predicates (orientation on the sphere and power test) reduce to orientation predicates on the data points in R 3 . Note that Section 3 shows that the weight of a point p is implicit, as it need not be explicitly computed throughout the entire algorithm.

Depending on the weights, some points can be hidden in a regular triangulation. We prove now that under some sampling conditions on the rounded data, there is actually no hidden point. Lemma 4.1. Let us assume that all data points lie within a distance δ from S. If the distance between any two points is larger than 2 √ Rδ, then no point is hidden. Proof. A point is hidden i it is contained inside the 3D convex hull of the set of data points S. Let p be a data point, at distance ρ from O. We have ρ ∈ [R -δ, R + δ]. Denote by d p the minimum distance between p and the other points. If

d p > (R + δ) 2 -ρ 2 , the set B(O, R + δ) \ B(p, d p )
is included in the half-space H + = {q : q -p, O -p > 0}. Under these conditions, all other points belong to H + and p is not inside the convex hull of the other points. It follows that if the distance between any two data points is larger than sup ρ (R + δ) 2 -ρ 2 = 2 √ Rδ, no point is hidden.

Let us now assume we use double precision oating point numbers as dened in the IEEE standard 754 [START_REF]IEEE standard for oating-point arithmetic[END_REF][START_REF] Goldberg | What every computer scientist should know about oating-point arithmetic[END_REF]. The mantissa is encoded using 52 bits. Let γ denote the worst error, for each Cartesian coordinate, done while rounding a point on S to the nearest point whose coordinates can be represented by double precision oating point numbers. Let us use the standard term ulp(x) denoting the gap between the two oating-point numbers closest to the real value x [START_REF] Muller | On the denition of ulp(x)[END_REF]. Assuming again that the center of S is O, one has γ ≤ ulp(R) = 2 -52+ log 2 (R) ≤ 2 -52 R. Then, δ in the previous lemma is such that δ ≤ 3/4γ < 2 -52 R. Using the result of the lemma, no point is hidden in the regular triangulation as soon as the distance between any two points is greater than 2 -25 R, which is highly probable in practice.

Implementation and Experiments

Both approaches presented in Section 4 have been implemented in C++, based on the cgal package that computes triangulations in R 2 . The package introduces an innite vertex in the triangulation to compactify R 2 . Thus the underlying combinatorial triangulation is a triangulation of the topological sphere. This allows us to reuse the whole combinatorial part of cgal 2D triangulations [START_REF] Pion | 2D triangulation data structure[END_REF] without any modication. However, the geometric embedding itself [START_REF] Yvinec | 2D triangulations[END_REF], bound to R 2 , must be modied by removing any reference to the innite vertex. A similar work was done to compute triangulations in the 3D at torus [START_REF] Caroli | Computing 3D periodic triangulations[END_REF][START_REF] Caroli | 3D periodic triangulations[END_REF], reusing the cgal 3D triangulation package [START_REF] Pion | 3D triangulation data structure[END_REF][START_REF] Pion | 3D triangulations[END_REF] as much as possible.

Also, the genericity oered in cgal by the mechanism of traits classes, that encapsulate the geometric predicates needed by the algorithms, allows us to easily use exactly the same algorithm with two dierent traits classes for our two approaches.

To display the triangulation and its dual, the code is interfaced with the cgal 3D spherical kernel [START_REF] Pedro | 3D spherical geometry kernel[END_REF][START_REF] Pedro | Design of the CGAL 3D Spherical Kernel and application to arrangements of circles on a sphere[END_REF], which provides primitives on circular arcs in 3D. The vertices of the triangulations shown are the projections on the sphere of the rounded-data points. The circular arcs are drawn on the surface of the sphere (see Figures 4, 6 and7). We compare the running time of our approaches with several available software packages on a MacBook Pro 3,1 equipped with a 2.6 GHz Intel Core 2 processor and 2GB 667 MHz DDR2 SDRAM7 (see Figure 5). We consider large sets of random data points8 (up to 2 23 points) on the sphere, rounded to double coordinates. Figure 6 indicates running times on some real-life data. Graph 1st of Figure 5 shows the results of our rst approach. We coded a traits class implementing the exact predicates presented in Section 4.1, together with semi-static and dynamic ltering [START_REF] Li | Recent progress in exact geometric computation[END_REF]. The non-linear behavior of the running time is due to the fact that our semi-static lters hardly ever fail for less than 2 13 points, and almost always fail for more than 2 18 points.

Graph 2nd shows the results of the second approach. One of the predened kernels9 of cgal provides us directly with an exact implementation of the predicates, ltered both semi-statically and dynamically. In our experiments we have observed that no point is hidden with such distributions, even when the data-set is large, which conrms in practice the discussion of Section 4.2.

The cgal 3D Delaunay triangulation (graph DT3) [START_REF] Pion | 3D triangulations[END_REF], with the same cgal kernel, also provides this convex hull as a by-product. We insert the center of the sphere to avoid penalizing this code with too many predicate calls on ve cospherical points that would always cause lters to fail.

For these three approaches, 3D spatial sorting reduces the running time of the location step of point insertion [START_REF] Delage | Spatial sorting[END_REF][START_REF] Buchin | Constructing Delaunay triangulations along spacelling curves[END_REF].

If the data points are lying exactly on a sphere, their Delaunay Triangulation can be extracted from an arrangement of geodesic arcs as computed by the code of Fogel and Setter [START_REF] Fogel | Software for Voronoi diagram on a sphere[END_REF][START_REF] Fogel | Exact implementation of arrangements of geodesic arcs on the sphere with applications[END_REF]. Since it is not the main purpose of their algorithm, the running times are not comparable: close to 600 seconds for 2 12 points. Note however that the code is preliminary and has not been fully optimized yet. No graph is shown.

We consider the following two software packages computing a convex hull in 3D,10 for which the data points are rst rounded to points with integer coordinates. Predicates are evaluated exactly using single precision computations. Graph HULL corresponds to the code [hul] of Clarkson, who uses a randomized incremental construction [START_REF] Clarkson | Four results on randomized incremental constructions[END_REF] with an exact arithmetic on integers [START_REF] Clarkson | Safe and eective determinant evaluation[END_REF].

Graph SUG shows the running times of Sugihara's code in Fortran [START_REF]Three-dimensional convex hulls[END_REF][START_REF] Sugihara | Laguerre Voronoi diagram on the sphere[END_REF].

Graph QHULL shows the performance of the famous Qhull package of Barber et al. [qhu] when computing the 3D convex hull of the points. The option we use handles round-o errors from oating point arithmetic by merging facets of the convex hull when necessary.

Renka computes the triangulation with an algorithm similar to our rst approach, but his software STRIPACK, in Fortran, uses approximate computations in double [START_REF] Renka | Algorithm 772: STRIPACK: Delaunay triangulation and Voronoi diagram on the surface of a sphere[END_REF]. Consequently, it performs quite well on random points (better than our implementations for small random data-sets), but it fails on some data-sets: Using STRIPACK, we did not manage to compute a triangulation of more than 2 19 random points (it returns an error ag). The same occurred for the inputs used to produce Figures 6 and7. Our implementations handle arbitrary data sets.

To test for exactness we devised a point set that is especially hard to triangulate because it yields many sliver-shaped triangles in the triangulation. This point set is dened as weather stations all around the world. Data and more information can be found at http://www.locationidentifiers.org/. Our second approach computes the result in 0.14 seconds, while Qhull needs 0.35 seconds, and the rst approach 0.57 seconds. STRIPACK fails on this data-set.

Conclusion

The results show that our second approach yields better timings than all the other tested software packages for large data-sets, while being fully robust. Note that it can in fact be used as well to compute the convex hull of points that are not close to a sphere: The center of the sphere can be chosen at any point inside a tetrahedron formed by any four non-coplanar data points.

The rst approach is slower but still one of the most scalable. It exactly computes the triangulation for input points with algebraic coordinates lying on the sphere, and thus ensures that in any case all points will appear in the triangulation. It is the only one to do so within reasonable time and thus being useful for real-world applications.

Figure 2 :

 2 Figure 2: Regular triangulation of a set of circles in the plane (their power diagram is shown dashed)

3p 1

 1 Remember that S is centered at O and has squared radius R 2 . = π S (c 1 )

Figure 4 :

 4 Figure 4: Delaunay triangulation (left) and Voronoi diagram (right) of the 9,031 post-oces in France (including Overseas Departments and Territories).

Figure 5 :

 5 Figure 5: Comparative benchmarks. The programs were aborted when their running time was above 10 minutes or in case of failure.

  show the Delaunay triangulation of S 250 and the Voronoi diagram of S 100 .

Figure 6 :

 6 Figure 6: Delaunay triangulation (left) and Voronoi diagram (right) of 20,950weather stations all around the world. Data and more information can be found at http://www.locationidentifiers.org/. Our second approach computes the result in 0.14 seconds, while Qhull needs 0.35 seconds, and the rst approach 0.57 seconds. STRIPACK fails on this data-set.

Figure 7 :

 7 Figure 7: Delaunay triangulation of S 250 (left), Voronoi diagram of S 100 (right). STRIPACK fails for e.g. n = 1, 500.

  Triangulation de Delaunay robuste et ecace, pour des points sur la sphère ou proche d'elle

Résumé :

RR n°7004

On S, the interior (respectively exterior) of a circle c that is not a great circle of S corresponds to the interior (respectively exterior) of the half-cone in 3D, whose apex is the center of S and that intersects S along c.INRIA

Further details: MAC OS X version 10.5.7; compiler g++ 4.3.2 with -O3 and -DNDEBUG ags, g77 3.4.3 with -O3 for Fortran. All running times mentioned exclude time used by input/output.

generated by CGAL::Random_points_on_sphere_3 INRIA

precisely CGAL::Exact_predicates_inexact_constructions_kernel RR n°7004

The plot does not show the results of the cgal 3D convex hull package[START_REF] Hert | 3D convex hulls[END_REF] because it is much slower than all other methods (roughly 500 times slower than Qhull).INRIA

http://www.inria.fr ISSN 0249-6399

Acknowledgments

We warmly acknowledge Ophir Setter and E Fogel who kindly worked on their code to give us access to it. We wish to thank Jean-Marc Schlenker for very interesting discussions on the de Sitter space.

This work was partially supported by the ANR (Agence Nationale de la Recherche) under the Triangles project of the Programme blanc ANR-07-BLAN-0319 http://www.inria.fr/geometrica/collaborations/triangles/.