Pascal Romon

Symmetries and conserved quantities for minimal surfaces

.

Jet bundles and differential equations

We follow the formalism in [START_REF] Olver | Applications of Lie groups to differential equations[END_REF] transposed in the fiber bundle setting. Though we are mostly interested in maps from a surface M to R 3 , we will write them as sections of the trivial bundle M × R 3 . This method follows the more general construction (see [START_REF] Saunders | The geometry of jet bundles[END_REF]) and explains the type of vector fields and transformations we will encounter.

Jet bundles and prolongation

Let N be a trivial bundle with base M m and fiber F n . Coordinates on M and F will be noted respectively x = (x i ) 1≤i≤m and u = (u α ) 1≤α≤n (called the independent and dependent variables). Two local sections s 1 , s 2 around p ∈ M are equivalent at p if they coincide up to first order, i.e. ∂ (u α •s 1 )

∂x i (p) = ∂ (u α •s 2 )
∂x i (p) for all i, α. The equivalence class is called the 1-jet of s 1 at p; the collection of all 1-jets is a manifold N 1 , called the first jet bundle, endowed with two bundle maps π 1 : N 1 -→ M and π 1,0 : N 1 -→ N . It is easy to describe N 1 using the coordinate u (1) = (x i , u α , u α i ) with u α i = ∂ (u α •s) ∂x i (p); projections π 1 and π 1,0 are obvious. One defines then a prolongation operator from Γ(N ) to Γ(N 1 ) by

pr 1 s(x) = (x i , s α , ∂ s α ∂x i ) i,α
In the same way, we can define the k-jet bundle N k -→ N -→ M and the k-th prolongation pr k : Γ(N ) -→ Γ(N k ). We need adopt a multi-index notation1 : I = (i 1 , . . . , i ) will denote a -tuple with 1 ≤ i 1 ≤ • • • ≤ i ≤ m and

∂ I = ∂ ∂x i 1 • • • ∂ ∂x i
Furthermore |I| = is called the length of I; if |I| = 0 then I = ∅ and ∂ I f = f . The k-jet at p of a local section is its equivalence class for the relation s 1 ∼ p s 2 if ∂ I s 1 (p) = ∂ I s 2 (p) for all |I| ≤ k. A coordinate on N k is u (k) = (x i , u α I ) for |I| ≤ k and any α, and the k-th prolongation is

pr k s(x) = (x i , ∂ I s α ) i,α,|I|≤k
We have a sequence of fiber bundles

• • • - → N k -→ • • • -→ N 1 -→ N ↓ ↓ ↓ • • • id -→ M id -→ • • • id -→ M id -→ M
and we define N ∞ -→ M to be its inverse limit. The infinite jet at p of a local section is its equivalence class for the relation s 1 ∼ p s 2 if ∂ I s 1 (p) = ∂ I s 2 (p) for all I. N ∞ is an infinite dimensional bundle with coordinate u (∞) = (x i , u α I ) α,I . The infinite prolongation pr ∞ is defined in a similar way.

Vector fields on N ∞ are sections of the tangent bundle T N ∞ ; they are written as linear combinations of the ∂ ∂x i and ∂ ∂u α I without limitation on the number of terms. Vector fields in X(N ) = Γ(T N ) can also be prolonged like the sections. However the formula is more delicate to establish and write; it involves the total derivatives D i ∈ Γ(T N ∞ ), defined in coordinates by

D i = ∂ ∂x i + α,I
≤ k). Then if X = i X i ∂ ∂x i + α U α ∂ ∂u α , pr k X = i X i ∂ ∂x i + α 0≤|I|≤k U α,I ∂ ∂u α I (k ∈ N ∪ ∞)
where

U α,I = i X i u α I∪i + D I U α - i X i u α i
One frequently uses a specific function algebra on N ∞ , namely F = lim

-→ F(N k ) = {f ∈ F(N ∞ ); ∃k, f ∈ F(N k )}.

Differential equations

A differential equation of order k in this context is just a closed subbundle R in N k ; locally it is given by the equations R(u (k) ) = 0 where R = (R 1 , . . . , R p ) is a p-tuple of functions in F(N k ). A solution is a section s ∈ Γ(E) such that pr k s ∈ R. We can also embed the problem in some superior jet bundle or even N ∞ ; we need then prolong R in the subbundle R ∞ locally defined by the set of equations

D K R j (u (∞) ) = 0 ∀K
In particular a variational principle is given by a lagrangian density2 L ∈ F(N k ). For any s ∈ Γ(N k ) set

L(s) = M L(pr k s(x)) dx 1 ∧ • • • ∧ dx m
The critical points of L are the sections s ∈ Γ(N ) whose prolongation prs stay in the subbundle determined by the Euler-Lagrange equations

E α (L)(u) = 0 for 1 ≤ α ≤ n (E α (L) generally belongs in F(N 2k )): E α (L) = I (-D) I ∂ L ∂u α I where (-D) I = (-1) |I| D I .
There eventually remains the problem of recognizing the immersions; that is not a priori implied by the previous study. Either one accepts solutions with singularities, or one restricts to the sections of the open set in N k given by: rank (u 1 , . . . , u m ) = m.

Symmetries and variational symmetries

Let G be a Lie group acting on the bundle N . For instance, G may be a group of transformations for the base space M , or fiber transformations; typically (x, u) → (x(x), ũ(u)). These are called classical transformations [START_REF] Bryant | Toward a geometry of differential equations[END_REF]. Elements of G may also be gauge transformations (x, u) → (x(x), ũ(x, u)); or even more generally point transformations (x, u) → (x(x, u), ũ(x, u)). Given a group G one computes its k-th prolongation, pr k G which acts on N k (k possibly infinite). Of course G may be only a local group of transformations. We may also consider the infinitesimal action of G, i.e. vector fields in the Lie algebra g, and we will always write g as a subalgebra of Γ(T N ).

Let R = (R 1 , . . . , R p ) be a partial differential equation of order k, i.e. a subbundle R in N k , satisfying some technical hypotheses 3 (clear for minimal surfaces). A group G leaves R invariant iff prG maps R into itself. A vector field X ∈ Γ(T N ) is an infinitesimal symmetry of R if e t prX maps R into itself, that is (pr k X)R = 0 on R (taking the Lie derivative). To check for such a symmetry we need only write the k-th prolongation and apply it to the functions R 1 , . . . , R p .

If R is the Euler-Lagrange equation for a Lagrangian L = M L, we define a variational symmetry to be a vector field

X = i X i ∂ ∂x j + α U α ∂ ∂u α such that (prX)L + LDiv (X 1 , . . . , X m ) = 0
for all u, where Div (X 1 , . . . , X m ) = D 1 X 1 +• • •+D m X m is called the total divergence of (X 1 , . . . , X m ). Such an infinitesimal symmetry preserves critical points of L hence (prX)E α (L) = 0 for all α. However this condition is too restrictive. Indeed it is enough for X to be a divergence symmetry : ∃P = (P 1 , . . . , P m ), (prX)L + LDiv (X 1 , . . . , X m ) = Div (P 1 , . . . , P m )

This ensures the invariance of the Euler-Lagrange4 equations under X.

Generalizations:

It is possible to look for more general type of infinitesimal symmetries and indeed we will meet some of these in the following. A first approach (due to Noether herself) consists in taking fields in F ⊗ Γ(T N ) :

X = i X i ∂ ∂x i + α U α ∂ ∂u α
with X i and U α in F (more generally X is a section of the pull back bundle π * ,0 (T N ) for some finite , where π ,0 : N -→ N ; see [START_REF] Saunders | The geometry of jet bundles[END_REF] for details). Another approach uses contact transformations and vector fields: we take X in Γ(T N ) instead of Γ(T N ); however we have to restrict to those which map solutions to solutions; there is then a compatibility condition written in terms of connection or differential ideal (see [START_REF] Bryant | Toward a geometry of differential equations[END_REF][START_REF] Tsujishita | On variation bicomplexes associated to differential equations[END_REF]). This might however take us too far from the classical transformations. Finally we may consider vector fields X defined only on R; that enlarges our scope (e.g. the conjugate family for minimal surfaces).

Noether's theorem

A conservation law is a (m -1)-form ω in Ω m-1 (M ) ⊗ F, whose restriction to R ∞ is closed for the total exterior derivative D. In other words, let s ∈ Γ(N ) be a solution with prolongation prs, then (Dω)(prs) = d(ω(prs)) = 0 Let γ be a (m -1)-chain in M ; then γ ω(prs) only depends on the homology class of γ, for any fixed solution s. We call this integral a conserved quantity.

Noether's theorem states the equivalence between variational symmetries and conservation laws. We will use it in the simple case where L is an order 1 lagrangian, and the variational symmetry is vertical, i.e. only changes the fiber N : X = X α ∂ ∂u α . Then using integration by parts:

(pr 1 X)L = α X α ∂ L ∂u α + i D i X α ∂ L ∂u α i = α X α ∂ L ∂u α - i X α D i ∂ L ∂u α i + i D i α X α ∂ L ∂u α i = XE(L) + i D i C i with C i = α X α ∂ L ∂u α i
We define the conservation law

ω = i (-1) i+1 C i dx 1 ∧ • • • ∧ dx i ∧ • • • ∧ dx m so that Dω = ( i D i C i ) dx 1 ∧ • • • ∧ dx m
which vanishes when evaluated at a prolonged solution u (∞) = prs. If X is a divergence symmetry, so that (pr 1 X)L = Div (P 1 , . . . , P m ), we need only replace C i by C i -P i .

Minimal surfaces

Let M be a simply-connected bounded domain of R 2 and N be the trivial bundle M × R 3 . We will see later how and when symmetries and conserved quantities extend to non-simply connected surfaces. For the sake of simplicity, we assume that M is the domain of a coordinate chart (x 1 , x 2 ). Minimal surfaces M -→ R 3 are the critical points of the order 1 lagrangian

L = M L( ∂ f ∂x 1 , ∂ f ∂x 2 ) dx 1 dx 2 with L ∈ F(N 1 ) is the map : u (1) → u 1 × u 2 ; × is the cross product in R 3
and . the Euclidean norm associated to the scalar product •|• . As the target space is a vector space, we use simpler index-free notations, like ∇ u I L for the vector

( ∂ L ∂u α I ) 1≤α≤3 .
Then the Euler-Lagrange equation in the second order jet bundle is:

E α (L)(u (2) ) = J (-D) J ∂ L ∂u α J (u (2) ) = - β u β 11 ∂ ∂u β 1 + u β 12 ∂ ∂u β 2 ∂ L ∂u α 1 + u β 12 ∂ ∂u β 1 + u β 22 ∂ ∂u β 2 ∂ L ∂u α 2
or in vector notation, calling ν the unit normal vector field, ν = (u 1 × u 2 )/L :

E(L)(u (2) ) = J (-D) J ∇ u J L(u (2) ) = -D 1 (u 2 × ν) -D 2 (ν × u 1 ) = -u 2 × D 1 ν -u 12 × ν -ν × u 12 -D 2 ν × u 1 = -u 2 × D 1 ν -D 2 ν × u 1 One checks that E(L)|u 1 = E(L)|u 2 = 0 thus E(L) = E(L)|ν ν = -2LHν = -2H(u 1 × u 2
) with H the mean curvature:

H = ν| u 2 2 u 11 + u 1 2 u 22 -2 u 1 |u 2 u 12 2 u 1 × u 2 2
If the immersion is conformal then H = 1 2L ∆u|ν . Furthermore ∆u is parallel to ν, hence E(L) = -∆u.

Symmetries

Notice that L does not depend on the x variable, but actually we know more : a change in parametrisation does not change the minimality 5 . This can be easily verified

: let X = X 1 ∂ ∂x 1 + X 2 ∂ ∂x 2 ; its first prolongation is (pr 1 X)L(u) = - α [(X 1 1 u α 1 + X 2 1 u α 2 ) ∂ L ∂u α 1 (u) + (X 1 2 u α 1 + X 2 2 u α 2 ) ∂ L ∂u α 2 (u)] = -X 1 1 u 1 + X 2 1 u 2 |∇ u 1 L -X 1 2 u 1 + X 2 2 u 2 |∇ u 2 L = -X 1 1 u 1 + X 2 1 u 2 |u 2 × ν -X 1 2 u 1 + X 2 2 u 2 |ν × u 1 = -(X 1 1 + X 2 2 ) ν|u 1 × u 2 = -LDiv (X 1 , X 2 )
This equation tells us that a change of variable x = φ(x) is a symmetry; indeed L may change, but the area L does not.

From now on we restrict our attention to infinitesimal symmetries keeping x fixed (so-called vertical generalized vector fields) 6 , so the term LDiv (X 1 , . . . , X m ) will always vanish, so we write:

X = 3 α=1 X α ∂ ∂u α
and X α ∈ F. The first prolongation is simply

pr 1 X = X + α [(D 1 X α ) ∂ ∂u α 1 + (D 2 X α ) ∂ ∂u α 2 ]
and noting

D i X = (D i X 1 , D i X 2 , D i X 3 ), X is a variational symmetry if (pr 1 X)L(u (1) ) = D 1 X|u 2 × ν + D 2 X|ν × u 1 = 0 (1)
Let us analyse what are the known symmetries.

Classical transformations: Classical vector fields7 leaving the metric in R 3

invariant will not change L, let alone H. These are simply translations X = constant (hence D 1 X = D 2 X = 0) and rotations in R 3 for which X(u) = A × u for some axis vector A.

(pr 1 X)L(u (1) ) = A × u 1 |u 2 × ν + A × u 2 |ν × u 1 = ν|(A × u 1 ) × u 2 + u 1 × (A × u 2 )
Jacobi's equality for the cross product implies:

(A × u 1 ) × u 2 + u 1 × (A × u 2 ) = (A × u 1 ) × u 2 + (u 2 × A) × u 1 = -(u 1 × u 2 ) × A hence (pr 1 X)L(u (1) ) = ν|A × (u 1 × u 2 ) = 0.
Dilations are not variational symmetries, however it will be useful to compute (pr 1 X)L(u (1) ) for X = u α ∂ ∂u α :

(pr 1 X)L(u (1) ) = u 1 |u 2 × ν + u 2 |ν × u 1 = Div (P 1 , P 2 )
with P 1 = u|u 2 × ν and P 2 = u|ν × u 1 ; indeed

Div (P 1 , P 2 ) = u 1 |u 2 × ν + u|u 12 × ν + u 2 × D 1 ν + u 2 |ν × u 1 + u|ν × u 12 + D 2 ν × u 1 = (pr 1 X)L(u (1) ) + u|u 2 × D 1 ν + D 2 ν × u 1 = (pr 1 X)L(u (1) ) -u|E(L)(u (1) )
which vanishes for minimal surfaces only.

Proposition 1 The classical variational vertical symmetries correspond exactly to the translations and affine rotations of R 3 .

Proof:

Let X ∈ X(R 3 ) be a variational symmetry; then (1) is true for any prolonged solution u = prs; in particular for conformal immersions, satisfying u 1 |u 2 = 0 and u 1 = u 2 . We get the simpler equation:

(pr 1 X)L(u (1) ) = dX(u) u 1 |u 1 + dX(u) u 2 |u 2 (2)
For any rotation R, Rs still is minimal and conformal. Fix any point p ∈ M , and choose an affine rotation keeping s(p) fixed (so u α (p) does not change); only u 1 and u 2 vary, and take actually arbitrary values (provided u 1 |u 2 = 0 and

u 1 = u 2 ). Clearly for any u ∈ R 3 , dX(u) is a linear map satisfying ∀ξ, η ∈ R 3 such that ξ|η = 0, ξ = η , dX(u) ξ|ξ + dX(u) η|η = 0 (3)
Then dX(u) is skew-symmetric (see lemma 2 below). That is not possible, except for constant X. Indeed the second order differential d 2 X(u) = d(dX)(u) is a symmetric vector valued form and we have both properties:

d 2 X(u) (ζ, η)|ξ = d 2 X(u) (η, ζ)|ξ symmetry d 2 X(u) (ζ, η)|ξ = d 2 X(u) (ζ, ξ)|η skew-symmetry of d 2 X(u)ζ
We get a contradiction by writing

d 2 X(u) (ζ, η)|ξ = -d 2 X(u) (ζ, ξ)|η = -d 2 X(u) (ξ, ζ)|η
so that the expression is both symmetric and skew-symmetric in ζ and η ! We conclude that d 2 X vanishes identically, so X(u) = dX(0) u+constant with dX(0) ∈ so(3), which is indeed an affine rotation (or simply a translation if dX(0) = 0).

Lemma 2 Let A ∈ L(R 3 ) be such that ∀ξ, η ∈ R 3 such that ξ|η = 0, ξ = η , Aξ|ξ + Aη|η = 0 then A is skew-symmetric.
Proof: Obviously skew-symmetry is a sufficient condition for that property to hold. Now let ξ be a unit real eigenvector of A with associated eigenvalue λ (there has to be at least one such). Then for any unit vector

η ∈ ξ ⊥ , Aη|η = -λ = -λ η|η . Thus on ξ ⊥ , ζ → Aζ|ζ is just the quadratic form -λ ζ 2 . Let (η 1 , η 2
) be a orthonormal basis of ξ ⊥ ; then applying the hypothesis yields:

0 = Aη 1 |η 1 + Aη 2 |η 2 = -2λ
So in the orthonormal basis (ξ, η 1 , η 2 ) the matrix of the bilinear form

A • |• is written   0 a b 0 0 c 0 -c 0   But reapplying the hypothesis to (η 2 , ζ+η 1 √ 2 ) yields: 0 = Aη 2 |η 2 + 1 2 A(ξ + η 1 )|(ξ + η 1 ) = 1 2 Aη 1 |ξ = a
The same computation gives b = 0. We conclude that A is skew-symmetric.

Let us now describe generalized symmetries.

The conjugate family:

We define an action of S 1 on the solutions by θ → u θ = cos θ u + sin θ v, where v is "the" conjugate immersion (since v is only defined up to translation, assume a point p ∈ M is fixed and and u(p) = v(p) = 0). If u is harmonic, v is just a conjugate harmonic map; otherwise there exists a change of coordinates φ such that û = u•φ is harmonic; then set v = v • φ -1 . It is straightforward to check that the definition does not depend on the choice of φ. Thus we define: X = d dθ |θ=0 u θ = v. A (somewhat long) calculation, eliminating φ, shows that

D 1 X = u 1 × ν D 2 X = u 2 × ν which corresponds to D 1 X = v 1 = -u 2 and D 2 X = v 2 = u 1 in the conformal case. Clearly condition (1) is satisfied. Although D i X exists for any u ∈ T N ∞ it is locally integrable if and only if D 2 D 1 X = D 1 D 2 X but D 1 (u 2 × ν) - D 2 (u 1 × ν) = -E(L)
vanishes exactly on minimal maps. So this is truly a generalized variational symmetry tangent only to the solution subbundle. Moreover, X =

x u 1 × ν dx 1 + u 2 × ν dx 2 cannot truly be described as a function on any k-th order jet space; but -considering analyticity -as a function of all derivatives, i.e. X ∈ F(N ∞ ) ⊗ X(N ).

The imaginary rotations: Start with the Weierstrass (hence conformal) representation for minimal surfaces, so u(x) = Re

x p (g -1 -g, i(g -1 + g), 2)η; following Lopez and Ros ( [START_REF] Pérez | Some uniqueness and nonexistence theorems for embedded minimal surfaces[END_REF]), we define an action of R * + : (g, η) → (λg, η), so

u λ (x) = Re x p (λ -1 g -1 -λg, i(λ -1 g -1 + λg), 2)η Thus X = du λ dλ |λ=1 = Re x -(g -1 + g), i(-g -1 + g), 0 η
Comparing it with the conjugate surface parametrisation v = Re

x p i(g -1 - g), -(g -1 + g), 2i η we infer that X = (v 2 , -v 1 , 0) = v × e 3 .
Obviously we can generalize to the field X = A × v for some fixed vector A (here again X is only defined on the solution subbundle). We get (pr 1 X)L(u (1)

) = A × (u 1 × ν)|u 2 × ν + A × (u 2 × ν)|ν × u 1 = A|(u 1 × ν) × (u 2 × ν) + (u 2 × ν) × (ν × u 1 ) = 2 u 1 × u 2 |A
Hence (pr 1 X)L(u (1) ) = Div (P 1 , P 2 ) with P 1 (u (1) ) = u × u 2 |A and P 2 (u (1) ) = u 1 × u|A . We see that its form is very similar to rotations (but for taking v instead of u); meaning of the name "imaginary rotations" will become more obvious later.

3 Applications of Noether's theorem

Formulation

For minimal surfaces it is very easy to define a Noether 1-form ω. Let X be a variational (possibly divergence) symmetry: (pr

1 X)L = D 1 P 1 + D 2 P 2 . Then set ω X = ( X|u 1 × ν + P 2 ) dx 1 + ( X|u 2 × ν -P 1 ) dx 2
and by construction ω X is closed, hence locally integrable on the solution subbundle into a function f X .

Translations and flux: Take X any constant vector A ∈ R 3 generating a translation subgroup, then there is no divergence term and

ω = A|(u 1 × ν)dx 1 + (u 2 × ν)dx 2
Clearly one can put together these results for various A by defining a

R 3 -valued 1-form ω = (u 1 × ν)dx 1 + (u 2 × ν)dx 2
then γ ω is a homological invariant, called the flux or force F along the curve γ. We note that the primitive of ω is nothing else than the conjugate map v (up to a constant); indeed the flux of u along γ is the period of v along that same curve.

Rotations and torque: For rotation with vector axis A, we get

ω = A × u|(u 1 × ν)dx 1 + (u 2 × ν)dx 2
Here again we define a vector-valued 1-form

ω = u × (u 1 × ν) dx 1 + u × (u 2 × ν) dx 2 .
The integral of ω along a curve is called the torque T of γ. It must be noted that unlike the flux, T (γ) changes when the immersion u is translated (because of the Lie bracket of the translation with the rotation); so if u -→ u+c then T (γ) -→ T (γ) + c × F(γ). A neat result can be proved using the torque:

Proposition 3 : Let M be a minimal surface with two (or three) vertical catenoidal ends, with axes labelled D 1 , D 2 (D 3 ). If there are two ends then both axes are identical. In there are three ends, then all three axes lie on the same vertical plane.

Proof: Let us first compute the flux and torque around a vertical catenoidal end 8 centered around the third coordinate axis. The end can be parametrized as an analytic graph :

x 3 (x 1 , x 2 ) = c log R + b + O(r -2 )
where c, b are constants and r = x = (x 1 ) 2 + (x 2 ) 2 . Then up to an additive factor O(r -3 )

u 1 =   1 0 cx 1 /r 2   , u 2 =   0 1 cx 2 /r 2   , ν =   -cx 1 /r 2 -cx 2 /r 2 1 -c 2 /2r 2  
8 meaning by that along a positively oriented curve representing the homology of the end.

u 1 ×ν =   c 2 x 1 x 2 /r 4 -1 + c 2 ((x 2 ) 2 -(x 1 ) 2 )/2r 4 -cx 2 /r 2   , u 2 ×ν =   1 + c 2 ((x 2 ) 2 -(x 1 ) 2 )/2r 4 -c 2 x 1 x 2 /r 4 cx 1 /r 2   so x =r ω = 2πc e 3 + O(r -2 )
Since the integral is independent of r for homological reasons, F = 2πc e 3 . For the torque note that it is hence invariant by vertical translation; we may as well assume b = 0. Then up to O(r -2 )

u × (u 1 × ν) =   c log r -c 3 log r((x 2 ) 2 -(x 1 ) 2 )/2r 4 -c(x 2 ) 2 /r 2 c 3 log r x 1 x 2 /r 4 + cx 1 x 2 /r 2 -x 1 + c 2 ((x 2 ) 2 -(x 1 ) 2 )x 1 /2r 4 -c 2 x 1 (x 2 ) 2 /r 4  
Hence the torque vanishes (taking the limit as r -→ +∞). Now if the catenoidal end is not centered on the x 3 -axis, but on an axis D, such that D -c is the x 3axis, then the torque around the end is T = c×F. If we have many catenoidal ends with respective fluxes F i = f i e 3 and torques

T i = c i × F i = f i c i × e 3
(taking all c i to be horizontal vectors). We choose oriented representatives γ i of all ends, and they bound a compact piece of the surface, hence the conserved quantity vanish. For instance, the sum of fluxes and the sum of torques vanish. For two ends, f 1 = -f 2 = 0 and

0 = T 1 + T 2 = f 1 (c 1 -c 2 ) × e 3
thus c 1 = c 2 . For three ends, the sum of torques is

3 1 f i c i × e 3 = 0
However we may translate the surface so that c 3 = 0; we then get f 1 c 1 +f 2 c 2 = 0, a non trivial linear relation. Of course similar relationship hold for a greater number of ends, though without giving sufficient information. This kind of relation is called a balancing formula.

Dilations: For dilations the conserved quantity may be computed (though dilations are no generalized variational symmetries), but it is simply a constant, indeed if X(u) = u, ω X = 0.

Conjugate family:

We denote here v a conjugate surface; notice however that v is not well-defined, since periods may be non-zero; so that v is defined on some cover of M ; furthermore v is only defined up to translation. For any closed curve γ, the invariant associated with the infinitesimal conjugate action is

C(γ) = γ v|dv A different choice of origin for v changes this invariant; if v - → v + c, C(γ) - → C(γ) + c|F(γ)
. which makes the invariant geometrically meaningful only for curves with zero flux 9 . Notice however that this integral is trivial for zero flux, yielding γ d v 2 = 0.

Lopez-Ros transformation: The infinitesimal symmetry is X = A × v for a chosen vector A; taking all direction in account we define here a vector valued 1-form ω by

ω = (v × v 1 -u × u 1 )dx 1 + (v × v 2 -u × u 2 )dx 2 = v × dv -u × du and we set P(γ) = γ ω. Clearly P is independent of translations of u, but if v - → v + c, then P - → P + c × F(γ).

Flat ends

We study in details a flat embedded annular end M . Assume for simplicity that M is horizontal. Let u(z) = Re z (1 -g 2 , i(1 + g 2 ), 2g)f dz be the standard Weierstrass representation, defined on a open neighbourhood of zero in C. Up to a change of orientation, we may assume g(0) = 0 with multiplicity k ≥ 2.

Embeddedness forces f to have a double pole. We may also choose the local coordinate around zero and rotate the surface so that g(z) = az k + O(z k+1 ) for some real number a > 0 and f (z) = e iα z -2 for some other real number α.

Then up to translation and noting z = re iθ

u(z) =   -r -1 cos(α -θ) + O(r 2k-1 ) r -1 sin(α -θ) + O(r 2k-1 ) 2a k-1 cos(α + (k -1)θ)r k-1 + O(r k )  
and the conjugate map is well-defined

v(z) =   -r -1 sin(α -θ) + O(r 2k-1 ) -r -1 cos(α -θ) + O(r 2k-1 ) 2a k-1 sin(α + (k -1)θ)r k-1 + O(r k )  
hence there is no flux around that end. Notice that the integration constants are chosen so that both ends are "centered". To imagine the end, we study its intersections with the asymptotic plane x 3 = 0, and that gives us θ ≡ ( π 2 -α)/(k-1)+O(r) mod π/(k-1) . So the intersection stays within bounded distance from 2(k -1) rays from the origin given by r → (-r -1 cos(αθ), r -1 sin(α -θ), 0) for any θ satisfying the previous relation. In particular for multiplicity 2 ends, the intersection stays close to the line R(sin 2α, cos 2α, 0).

The positively oriented curve γ will be for constant r θ → re -iθ (coherent with the Gauss map pointing towards the south pole). Let us compute the torque T = γ u × dv.

u× ∂v ∂θ =   -r -1 cos(α + θ) + O(r 2k-1 ) r -1 sin(α + θ) + O(r 2k-1 ) 2a k-1 cos(α -(k -1)θ)r k-1 + O(r k )   ×   -r -1 cos(α + θ) + O(r 2k-1 ) r -1 sin(α + θ) + O(r 2k-1 ) -2ar k-1 cos(α -(k -1)θ) + O(r k )   =   -2ka k-1 r k-2 sin(α + θ) cos(α -(k -1)θ) + O(r k-1 ) -2ka k-1 r k-2 cos(α + θ) cos(α -(k -1)θ) + O(r k-1 ) O(r 2k-2 )  
When we integrate as r -→ 0, all positive powers of r vanish and only constant terms may remain, that is in r k-2 iff k = 2; otherwise the torque simply vanishes 10 . Let us compute T in the k = 2 case. Then

T = -4a π -π   sin(α + θ) cos(α -θ) cos(α + θ) cos(α -θ) 0   dθ = -4πa   sin 2α cos 2α 0  
That is a horizontal vector, depicting the asymptotic direction characterizing the end. Indeed the intersection with the asymptotic horizontal plane is parametrised by : r → ±r -1 (sin 2α, cos 2α) + O(1); hence there is a welldefined asymptotic direction, with vector T . The sign of T determines on which side of this asymptotic line is the surface above (or under) the asymptotic plane (see figure ??).

[[ Image ]]

That stays true, even under less restrictive conditions as those used above for the computations. In the general case, the end is still horizontal but we do not assume that the leading coefficients in g and f expansion are respectively real and unitary. Then

T = 4π   Im c Re c 0   where c = Res ( √ g f dz; 0) 2
Indeed it is possible to choose locally a square root h of g defined up to sign. So c is a well defined invariant.

The conserved quantity for the infinitesimal conjugation is

C = γ v|dv v| ∂v ∂θ =   -r -1 sin(α + θ) + O(r 2k-1 ) -r -1 cos(α + θ) + O(r 2k-1 ) 2a k-1 sin(α -(k -1)θ)r k-1 + O(r k )   |   -r -1 cos(α + θ) + O(r 2k-1 ) r -1 sin(α + θ) + O(r 2k-1 ) -2ar k-1 cos(α + (1 -k)θ) + O(r k )   = O(r 2k-2 )
which vanishes when integrated along γ.

For the Lopez-Ros transformation, P = γ v × dv -u × du. Here again the only interesting terms are those in r 0 , requiring k = 2; let us assume that and single out the relevant terms

v × ∂v ∂θ =   -r -1 sin(α + θ) + O(r 3 ) -r -1 cos(α + θ) + O(r 3 ) 2ar sin(α -θ) + O(r 2 )   ×   -r -1 cos(α + θ) + O(r 3 ) r -1 sin(α + θ) + O(r 3 ) -2ar cos(α -θ) + O(r 2 )   = 2a   cos(α + θ) cos(α -θ) -sin(α + θ) sin(α -θ) -cos(α + θ) sin(α -θ) -sin(α + θ) cos(α -θ) 0   +
P = 8πa   cos 2α -sin 2α 0   = 2   -1 1 0   T
In the general case, let ν the unit normal at the end, then

P = 2ν × T Proposition 4
Let M be a singly-periodic complete minimal cylinder, with two flat ends in the quotient, and suppose M is symmetric with respect to a vertical plane P . Then both ends are asymptotically identical.

Proof: Let T be the period and γ a generator of the homology of the cylinder. Then both T and F(γ) belong in P the plane of symmetry. Let γ + T be the curve γ translated by T ; then T (γ + T ) = T (γ) + T × F(γ). So that if T i denotes the torque around the i-th end

T (γ + T ) -T (γ) = T × F(γ) = T 1 + T 2
And by symmetry, both torques T i are perpendicular to P . Consider now the Lopez-Ros invariant P. Since the conjugate map is not well defined (there may be flux along γ), P is not well-defined; however P = P|F(γ) is. Translation in u leaves P invariant, thus P (γ + T ) = P (γ). So we must have at the ends

0 = P 1 + P 2 = 2ν 1 × T 1 + 2ν 2 × T 2 = 2ν 1 × (T 1 -T 2 )
where ν i is the normal at then end, and belongs in P . Hence T 1 = T 2 . Both ends are indeed asymptotically equal.

If there is no plane of symmetry, we can however say something. Suppose T × F(γ) = 0 (can it vanish ?). Define then the plane P = Vect(T, F(γ)). This plane is vertical. Indeed T (γ + T ) -T (γ) = T × F(γ) = T 1 + T 2 which is horizontal. Use the splitting R 3 = P ⊕ P ⊥ to write T i = T P i ⊕ T ⊥ i . If F(γ) is not vertical we can follow the same reasoning as above, and conclude that T ⊥ 1 = T ⊥ 2 , and then T P 1 = T P 2 . What remains is P 1 + P 2 = 2ν 1 × T P 1 which is equal to ±F(γ) × F(δ) where δ is the other generator of the torus (in the quotient).

for any prolonged solution φ; and the condition for keeping conformality is simply Z|φ = 0. The simplest possibility is for linear and zero-order vector fields: Z(φ) = Aφ for constant matrix A. Then conformality imposes that A belong in CI 3 ⊕ so(3, C). These are the already known transformations: I 3 is the dilation, iI 3 is the conjugation, so(3, R) corresponds to rotations, and i so(3, R) to the 3-dimensional family of imaginary rotations (Lopez-Ros transformations).

How to get the corresponding real field:

Let Z be a infinitesimal symmetry. The corresponding variational symmetry Ẑ is simply Ẑ(u) = Re x x 0 Z(φ) dx; of course it is defined up to translation, and is by construction well defined locally (there may be period problems). The symmetry Ẑ may be a divergence symmetry, indeed

(pr Ẑ)L(u) = D 1 Ẑ|u 1 + D 2 Ẑ|u 2 = X|u 1 -Y |u 2 = D 1 X|u -D 1 X|u -D 2 Y |u + D 2 Y |u = D 1 P 1 + D 2 P 2 + u|D 2 Y -D 1 X = Div P with P 1 = X|u , P 2 = -Y |u , since holomorphicity implies D 2 Y = D 1 X.
On the other hand

(pr Ẑ)L(u) = Re [Z(φ)]|Re [φ] + Im [Z(φ)]|Im [φ] = Re Z(φ)| φ = 1 2 Z. φ| φ
Hence Ẑ is a (non-divergence) variational symmetry iff Z keeps the metric; that is obvious: indeed for conformal maps, fixing the metric is equivalent to fixing the area. If Z fixes the metric and depends only on φ, then it falls under the classification in [START_REF] Lawson | Lectures on minimal submanifolds[END_REF]. It is clear that a metric-preserving vector field in C [START_REF] Hélein | Applications harmoniques, lois de conservation et repères mobiles[END_REF] has to be linear (??). The dilation and the imaginary rotations however are truly divergence symmetries.

The Noether 1-form can be written as

ω = -( Ẑ|u 2 + Y |u )dx 1 + ( Ẑ|u 1 -X|u )dx 2
However in this definition there may integration constants. So the conserved quantity is not always well-defined.

If Ẑ - → Ẑ + c then ω - → ω -c|u 2 dx 1 + c|u 1 dx 2 so γ ω - → γ ω + c|F(γ) If u - → u + c then ω - → ω + D 2 Ẑ|c dx 1 -D 1 Ẑ|c dx 2 so γ ω - → γ ω + c| γ D 2 Ẑ dx 1 -D 1 Ẑ dx 2
We may write ω in terms of the complex variable z. Fixing an arbitrary integration starting point, and denoting W = Z and Φ

= φ ω = -Re Z| -Im φ + -Im Z|Re φ dz + dz 2 + Re Z|Re φ -Re Z|Re φ dz -dz 2i = 1 8i W + W |φ -φ + Z -Z|Φ + Φ (dz + dz) + 1 8i W + W |φ + φ -Z + Z|Φ + Φ (dz -dz) = 1 4i W + W |φ -Z|Φ + Φ dz + 1 4i -W + W | φ + Z|Φ + Φ dz So ω = 1 2 Im α with α = ( W + W |φ -Z|Φ + Φ )dz.
Since ω is a homological invariant, we may integrate along a curve |z| = r = c te in the z-plane, and the result has to be independent from r.

Examples of non linear fields: Z(φ) = hφ for any C-valued holomorphic map h (for instance h = φ α ). More interesting: for any non zero integer k the generalized field

Z = D k φ × φ = φ k.1 × φ. Clearly φ k.1 × φ|φ = 0. For instance k = 1 and (pr Ẑ)L = u 11 × u 1 -u 12 × u 2 |u 1 + u 11 × u 2 + u 12 × u 1 |u 2 = 2 u 12 |u 1 × u 2
which does not vanish in general, so area is not fixed by Ẑ. Let us analyse this case in detail. For instance for a catenoidal end at z = 0 let φ = (1g 2 , i(1 + g 2 ), 2g)f and

g(z) = z(a + bz + cz 2 + O(z 3 )) f (z) = z -2 (α + γz 2 + O(z 3 )) So g 2 = z 2 (a 2 + 2abz + (b 2 + 2ac)z 2 + O(z 3 )). Then φ = z -2   α + (γ -a 2 )z 2 + O(z 3 ) i(α + (γ + a 2 )z 2 + O(z 3 )) 2z(aα + αbz + (aγ + αc)z 2 + O(z 3 ))   φ =   -2αz -3 + O(1) -2iαz -3 + O(1) -2aαz -2 + 2(γ + αc) + O(z)   Φ =   -αz -1 + (γ -a 2 )z + O(z 2 ) -iαz -1 + i(γ + a 2 )z + O(z 2 ) 2aα log z + 2αbz + (aγ + αc)z 2 + O(z 3 )   Z = φ × φ = αz -2   2iaαz -2 + 4ibαz

Lie brackets (unfinished)

The Lie bracket between generalized vertical symmetries is defined by (see [START_REF] Olver | Applications of Lie groups to differential equations[END_REF]) where (α, λ) are defined up to sign. For any f : ∆ -→ R 3 there exists φ ∈ D such that f • φ is conformal; φ is unique up to left composition by C; if ones imposes the 3-dim following conditions φ(0) = 0 ∂ 1 φ 2 (0) = 0 then φ is unique. We call D + the group of such diffeomorphisms.

First order approximation

Let u(ε) = u + εv be a local deformation (dropping here and afterwards all higher order terms); there exists diffeomorphisms φ(ε) (unique if chosen in

D + ) such that u(ε) • φ(ε) is conformal. Let Y = d dε φ(ε) |ε=0 u(ε) • φ(ε) = (u + εv) • (φ + εY ) = u • φ + ε(v • φ + Du • φ.Y ) ∂ 1 (u(ε) • φ(ε)) = (u • φ) 1 + ε{(v • φ) 1 + (Du • φ) 1 .Y + Du • φ.Y 1 } ∂ 2 (u(ε) • φ(ε)) = (u • φ) 2 + ε{(v • φ) 2 + (Du • φ) 2 .Y + Du • φ.Y 2 } First conformality condition 1 2ε ( ∂ 1 (u(ε)•φ(ε)) 2 -∂ 2 (u(ε)•φ(ε)) 2 ) = (u•φ) 1 |(v•φ) 1 +(Du•φ) 1 .Y +Du•φ.
(v • φ + Du • φ.Y ) 1 2 -(v • φ + Du • φ.Y ) 2 2 = 0 (v • φ + Du • φ.Y ) 1 |(v • φ + Du • φ.Y ) 2 = 0 + directness !

  and D I = D i 1 • • • D i (when working in N k one simply writes a finite sum with |I|

r - 1

 1 1 cos(α + θ) + O(r 3 ) r -1 sin(α + θ) + O(r 3 ) 2ar cos(α -θ) + O(r 2 ) sin(α + θ) + O(r 3 ) r -1 cos(α + θ) + O(r 3 ) 2ar sin(α -θ) + O(r 2 ) + θ) sin(α -θ) -cos(α + θ) cos(α -θ) cos(α + θ) sin(α -θ) + sin(α + θ) cos(α -θ)

- 3 - 2 3 6 Élucubrations

 326 -1 + 2i(aγ + 3αc -a 3 + γ) + O(z) -2aαz -2 -4bαz -1 -2(aγ + 3αc + a 3 + γ) + O(z) 4ia 2 z -1 + O(1) 2ibαz -2 -2i(aγ + 3αc -a 3 + γ)z -1 + O(log z) aαz -3 + 2bαz -2 + 2(aγ + 3αc + a 3 + γ)z -1 + O(log z) -2ia 2 z -2 + O(z -1 )La définition d'Olver pour les symétries variationnelles (divergente ou non) réclame que l'équation soit vraie pour tout u, et pas seulement pour les prolongations des solutions. Ce n'est pas le cas par contre des symétries des EDP, qui doivent simplement annuler l'expression pour toute solution ! Il semble cependant que pour appliquer Noether, il suffise d'avoir une symétrie divergente restreinte aux solutions (soit prX(L) + LDiv ξ = Div P seulement sur les solutions d'Euler-Lagrange) !! Alors la dilatation serait bien génératrice d'une loi de conservation (quoique triviale).

[ 8 Finding the conformal parametrization 8 . 1

 881 [X, Y ]] = prX(Y ) -prY (X) = I α,β (D I X β ) , Y ]] = I (d u I Y D I X -d u I X D I Y ) For classical vector fields, [[X, Y ]] = [X, Y ].The space of diffeomorphismsLet ∆ be the unit disk and D (resp C) be the groups of orientation preserving (conformal) diffeomorphisms. Notice that C is a 3-dimensional Lie group consisting in maps z → λz + α ᾱz + λ |λ| = 1, |α| < 1

Y 1 -

 1 (u • φ) 2 |(v • φ) 2 + (Du • φ) 2 .Y + Du • φ.Y 2 = 0 Second conformality condition ε -1 ∂ 1 (u(ε)•φ(ε))|∂ 2 (u(ε)•φ(ε)) = (u•φ) 1 |(v •φ) 2 +(Du•φ) 2 .Y +Du•φ.Y 2 + (u • φ) 2 |(v • φ) 1 + (Du • φ) 1 .Y + Du • φ.Y 1 = 0 which tells us that the basis ((v •φ+Du•φ.Y ) 1 , (v •φ+Du•φ.Y ) 2) is (almost) conformal (and direct):

Warning: this notation (as found in[START_REF] Olver | Applications of Lie groups to differential equations[END_REF]) is standard in PDE theory but not in jet bundle theory (see[START_REF] Saunders | The geometry of jet bundles[END_REF]).

if N is not trivial, L is rather defined as a section of the pull-back bundle π * k ( m T * M ), where π k is the projection N k -→ M .

[START_REF] Hélein | Applications harmoniques, lois de conservation et repères mobiles[END_REF] we ask for local solvability, which is implied by analyticity and normality of the Euler-Lagrange equations.

The Euler-Lagrange equations usually have higher order than L; hence this condition is a simpler one; however there may be symmetries of the Euler-Lagrange equation which are not variational symmetries (even not divergence symmetries).

for that very reason we often think of the minimal surface independently of its parametrisation.

as proven in[START_REF] Olver | Applications of Lie groups to differential equations[END_REF], we can do so without loss of generality.

in the sense defined above: vector fields on R 3 only.

In this latter case only, the starting point in the integral does not matter.

A few simple consequences may be drawn from this and the previous proposition, such as: surfaces with two catenoidal ends and flat ends of multiplicities greater than 2 have the same axis for both catenoidal ends (e.g. Costa's surface).

Other symmetries

The harmonic+conformal subbundle:

For simplicity reasons, we may consider instead of the minimal surface equation, the harmonic+conformal equation, defined by

The corresponding subbundle H ⊂ N 2 is obtained through (4) and its derivatives, up to order 2. This adds in theory 4 equations

Using the condition u 11 + u 22 = 0, we see that two of the previous conditions are redondant, hence we get a rank 7 system

Hence H is a dimension 13 subbundle of N 2 .

Harmonic conformal analysis

Let us now consider only minimal immersions f which are conformal (hence harmonic). The equation reduces to two (non-variational) PDEs, namely the Cauchy-Riemann equations for

and the conformality condition G(φ) = φ|φ = α (φ α ) 2 = 0. Notice that by using this formulation one avoids the problem of choosing an integration constant; two solutions differing by a translation are identified. Also firstorder symmetries in u become zero-order in φ (as will be illustrated below). However it must be remembered that φ does not actually define the minimal immersion, only φ dz does.

Using complex notations as a shortcut, let Z = X + iY be an infinitesimal symmetry acting on the prolonged φ-space. Then Z preserves F iff Z is holomorphic in its variables:


Christophe Lecoutre
email: lecoutre@cril.fr

Chavalit Likitvivatanavong
email: chavalit@comp.nus.edu.

Scott G Shannon

Roland H C Yap
email: ryap@comp.nus.edu.

Yuanlin Zhang
email: yzhang@cs.ttu.edu

M R C Van Dongen

Maintaining Arc Consistency with Multiple Residues

Keywords: Arc Consistency, Residual Supports, MAC

Exploiting residual supports (or residues) has proved to be one of the most cost-effective approaches for Maintaining Arc Consistency during search (MAC). While MAC based on optimal AC algorithm may have better theoretical time complexity in some cases, in practice the overhead for maintaining required data structure during search outweighs the benefit, not to mention the more complicated implementation. Implementing MAC with residues, on the other hand, is trivial.

In this paper we extend previous work on residues and investigate the use of multiple residues during search. We first give a theoretical analysis of residue-based algorithms that explains their good practical performance. We then propose several heuristics on how to deal with multiple residues. Finally, our empirical study shows that with a proper and limited number of residues, many constraint checks can be saved. When the constraint check is expensive or a problem is hard, the multiple residues approach is competitive in both the number of constraint checks and cpu time.

Introduction

Maintaining Arc Consistency (MAC) [START_REF] Sabin | Contradicting conventional wisdom in constraint satisfaction[END_REF] has been considered one of the most efficient algorithm for solving large and hard constraint satisfaction problems. At its core is the Arc Consistency algorithm (AC), whose efficiency plays a vital role in the overall performance of MAC. [START_REF] Lecoutre | A study of residual supports in arc consistency[END_REF] and [START_REF] Likitvivatanavong | Arc consistency during seach[END_REF] show that when arc consistency is enforced during search (ACS), a MAC3-like algorithm that simply reuses supports found earlier, called residual supports or residues, can outperform the optimal algorithm MAC2001/3.1 [START_REF] Bessière | An optimal coarsegrained arc consistency algorithm[END_REF]. Since a single residue is good, this paper investigates the use of multiple residues. We give a theoretical explanation for the effectiveness of residues. We then investigate various heuristics for using multiple residues. Finally, we perform extensive experiments with different variations of multiple residues. Our results show that when the number of residues is small (1 to 5), the number of constraint checks decreases sharply with a moderate increase of extra cost like validity checks.

Preliminaries

A (finite) Constraint Network (CN) P is a pair (X , C ) where X is a finite set of n variables and C a finite set of e constraints. Each variable X ∈ X has an associated domain containing the set of values allowed for X. The initial domain X is denoted by dom init (X); the current one by dom(X). Each constraint C ∈ C involves an ordered subset of variables of X called scope (denoted by scp(C)), and an associated relation (denoted by rel(C)). For each r-ary constraint C with scp(C) = {X 1 , . . . , X r }, rel(C) ⊆ r i=1 dom init (X i ). For any t = (a 1 , . . . , a r ) of rel(C), called a tuple, t[X i ] denotes a i . With a total order on the domains, tuples can be ordered using a lexicographic order ≺. To simplify the presentation, we use two special values ⊥ and ⊤ where ⊥ ≺ t ≺ ⊤ for any tuple t.

Let C be an r-ary constraint and scp(C) = {X 1 , . . . , X r }, an r-tuple t of r i=1 dom init (X i ) is said to be: (1) allowed by C iff t ∈ rel(C), (2) valid iff ∀X i ∈ scp(C), t[X i ] ∈ dom(X i ), (3) a support in C iff it is allowed by C and valid, and (4) a conflict iff it is not allowed by C and valid. A tuple t is a support of (X i , a) in C when t is a support in C and t[X i ] = a. A constraint check determines if a tuple is allowed. A validity check determines if a tuple is valid. A solution to a constraint network is an assignment of values to all the variables such that all the constraints are satisfied.

A pair (X, a), with X ∈ X and a ∈ dom(X), is generalized arc-consistent (GAC) iff ∀C ∈ C where X ∈ scp(C), there exists a support of (X, a) in

C. P is GAC iff ∀X ∈ X , dom(X) = ∅ and ∀a ∈ dom(X), (X, a) is GAC. For binary constraint networks GAC is called AC. A CN-value is a triplet (C, X, a) where C ∈ C , X ∈ scp(C), a ∈ dom(X).

Residual Supports

We illustrate the concept of residual support with GAC3. A residue for a CN-value is a support that has been previously found and stored for future use. Unlike the last structure in AC2001/3.1, a residue for a value might not be a lower bound of the current supports of the value. The concept of residue has been introduced under its multi-directional form by [START_REF] Lecoutre | Exploiting multidirectionality in coarse-grained arc consistency algorithms[END_REF] and under its uni-directional form by [START_REF] Likitvivatanavong | Arc consistency in MAC: a new perspective[END_REF].

GAC3 with Residual Supports

GAC3 with residual supports is shown in Algorithm 1. The algorithm adds a three-dimensional array res, which is initialized to ⊥. For a CN-value (C,X,a), res[C, X, a] stores the residue for (X, a) with respect to C.

Algorithm 1: GAC3 rm (P = (X , C ) : Constraint Network) for each C ∈ C ∧ X ∈ scp(C) ∧ a ∈ dom(X) do res[C, X, a] ← ⊥ 1 Q ← {(C, X) | C ∈ C ∧ X ∈ scp(C)} 2 while Q = ∅ do 3 extract (C, X) from Q 4 if revise(C, X) then 5 Q ← Q ∪ {(C ′ , X ′ ) | C ′ ∈ C ∧ C ′ = C ∧ X ′ = X ∧ {X, X ′ } ⊆ scp(C ′ )} 6 4
In each revision of arc (C, X) by revise (Algorithm 2), the validity of the residue for each CNvalue (C, X, a) is tested first (line 3). If it fails, a new support is searched from scratch (line 4). If a support t is found, multi-directionality is exploited to update the residues of all values of t (line 8), since t is also a support of t[Y ] for all Y ∈ scp(C). Consequently, r -1 residues (r is the arity of C) are obtained for other values in the tuple with no effort. The uni-directional form, in contrast, only updates the residue of the CN-value (C, X, a).

Algorithm 2: revise(C: Constraint, X: Variable): Boolean , t[X]]. Function isValid (Algorithm 3) determines whether or not the given tuple is valid (⊥ is not valid). Function seekSupport (Algorithm 4) determines from scratch a support for (X, a) in C. It uses function setNextValid (not listed) which returns either the smallest valid tuple t ′ built from C such that t ≺ t ′ and t ′ [X] = a, or ⊤ if it does not exist. ⊤ does not belong to any relation.

nbElements ← |dom(X)| 1 for each a ∈ dom(X) do 2 if isValid(C,res[C,X,a]) then continue 3 t ← seekSupport(C, X, a) 4 if t = ⊤ then 5 remove a from dom(X) 6 else 7 for each Y ∈ vars(C) do addResidue(C, Y, t) 8 return nbElements = |dom(X)| 9 Algorithm 3: isValid(C: Constraint, t: Tuple): Boolean if t = ⊥ then return false 1 for each X ∈ scp(C) do 2 if t[X] / ∈ dom(X)

Complexity Results

To understand why residues work, we present some results for binary problems. In particular, we study the complexity of AC3 rm when used stand-alone and when embedded in MAC. Without any loss of generality, we assume that each domain contains exactly d values.

Proposition 1 AC3 rm has a worst-case space complexity of O(ed) and a worst-case time complexity of O(ed 3 ).

Proof

We assume here that each constraint is represented in an intentional form with space complexity O(1), and n < e (otherwise, each connected component of the network can be analyzed similarly). The space required for a CN is O(e + nd), the space for Q is O(e * 2) = O(e), and the space for res is O(e * 2 * d) = O(ed). Thus, space complexity of AC3 rm is O(ed). Like the optimality proof of AC3.1 [START_REF] Bessière | An optimal coarsegrained arc consistency algorithm[END_REF], the number of validity checks performed by AC3 rm will not exceed the number of constraint checks performed by AC3. As AC3 requires O(ed 3 ) in the worst case, so is AC3 rm .

The analysis can be refined to consider the tightness of the constraints ( |allowed tuples| |possible tuples| ).

Definition 2 A constraint C is tightness-bounded iff for any CN-value involving C, either its number of supports is O(1) or its number of conflicts is O(1) when d → ∞.

Many common constraints are tightness-bounded. For example, X = Y or X = Y . For equations, each value is supported at most once. For dis-equations, each value allows at most one conflict. In practice, we observe that AC3 rm behaves in an optimal way when applied to constraints of small or high tightness.

Proposition 3 Applied to a constraint network involving tightness-bounded constraints, AC3 rm admits a worst-case time complexity of O(ed 2 ), which is optimal.

Proof In [START_REF] Lecoutre | A study of residual supports in arc consistency[END_REF], it is shown that the worst-case accumulated time complexity of seekSupport for a CN-value (C, X, a) is O(cs + d) where c is the number of conflicts of (X, a) in C and s the number of supports of (X, a Proposition 3 shows that AC3 rm behaves optimally when constraints are tightness-bounded. This suggests AC3 rm should be quite competitive, compared to optimal algorithms like AC2001/3.1, on highly structured problems. These results are confirmed when the state-of-the-art generic algorithm MAC [START_REF] Sabin | Contradicting conventional wisdom in constraint satisfaction[END_REF] is considered. The following result is directly obtained from previous propositions and the fact that AC3 is an incremental algorithm so no maintenance of data structures is necessary when backtracking. MAC3 rm denotes MAC embedding AC3 rm .

Proposition 4 MAC3 rm admits a worst-case space complexity of O(ed), and for any branch of the search tree admits: (1) a worst-case time complexity of O(ed 3 ), and (2) a worst-case time complexity of O(ed 2 ) for a constraint network involving tightness-bounded constraints.

These theoretical results partially justify the data in [START_REF] Likitvivatanavong | Arc consistency in MAC: a new perspective[END_REF][START_REF] Lecoutre | A study of residual supports in arc consistency[END_REF]. Further, they hold for single as well as (a bounded number of) multiple residues.

We will focus on the practical aspects of multiple residues from now.

Fundamentals of Multiple Residues

Given the results of the single residue approach [START_REF] Lecoutre | A study of residual supports in arc consistency[END_REF][START_REF] Likitvivatanavong | Arc consistency in MAC: a new perspective[END_REF], it is interesting to see if we can achieve better performance by using multiple residues.

Using more than one residue allows us to record not only the latest support, but also a selected subset of the past supports. The memory space for holding these residues for each CN-value is denoted residue store.

To find a support, an AC or ACS algorithm with multiple residues always considers the residue store first. If all residues fail the validity test, the algorithm looks for a new support which may then be added to the store. When the store is full, ideally the residue that has the least chance of being a support in the future is deleted. We approximate the idealized case using a number of heuristics. They can be used to evaluate each residue and return a value called utility that approximates the idealized probability. The utility also guides in identifying which residues to check first when looking for a support.

We present a generic algorithm and several policies that define the utility of a residue.

A Generic Multiple Residue Algorithm

The generic multiple residue algorithm is similar to the single one. The difference lies in how the residue store is updated and how the validity of residues is checked. Function addResidue is redefined in Algorithm 5 and isValid (used in Algorithm 2) is replaced by existValid (Algorithm 6).

There are three algorithm variants differing in management of utilities. The static variant computes a utility score only when a new support is added to the residue store. The dynamic and fully dynamic variants update the utility of the first found valid residue (line 4 of Algorithm 6). For the fully dynamic variant, the score of a residue is updated whenever its validity is checked (line 6 of Algorithm 6). The different variants are differentiated by the global variable updateType.

Algorithm 5: addResidue(C: Constraint, X: Variable, t: Tuple) The following data structures implement the ideas. Residue stores have fixed size maxR. The utility function is denoted by f . res[C, X, a] is now a Dictionary S with the following operations. max(S) returns a pair (t, u) where the residue t has the highest utility score u; pred(S, (t, u)) returns a pair (t ′ , u ′ ) where t ′ is the residue of utility u ′ which is the next smaller score than u, and returns ⊥ when there is none; insert(S, (t, u))/delete(S, (t, u)) inserts/deletes the pair (t, u); update(S, (t, u ′ )) updates the tuple (t, u) in the store by replacing u with u ′ and updates the tuple's position in S. size(S) returns the number of residues in S. In practice, a Dictionary can be implemented using a variety of data structures, for example balanced trees. We remark that the new support might not automatically become a residue when the residue store is full. For a new support to become a residue, its utility score must improve on the minimum score of the store (line 4 of Algorithm 5).

(generic routine) u ← f (t) 1 if size(res[C, X, t[X]]) = maxR then 2 (t0, u0) ← min(res[C, X, t[X]]) 3 if u0 < u then delete(res[C, X, t[X]], (t0, u0))) 4 if size(res[C, X, t[X]]) < maxR then insert(res[C, X,

Policies Defining the Utility Function

We consider the following policies which are based on evaluating the probability (utility) of a support being valid in the future:

1. Level of the search tree. The rationale is that the support found at a deeper level should be more robust to change in the network caused by a new assignment or backtracking. When the search backtracks, previously removed values are restored, but the current support is still valid until the search tree branches off to a different path. LEVELMIN policy defines the utility of a new support or a residue in the store as the current level of the search tree when it was found or checked (for validity).

2. Domain size. The rationale is that the support found when the domain is smaller is more robust. DOMMIN policy defines the utility of a support (or residue) based on the size of the relevant domains when it was found (or its validity is checked). Specifically, if revise(C, X) is invoked and the tuple t is found as a new support of a ∈ X, then

f (t) = Y ∈scp(C)-X |dom(Y )|.
3. Chronology. The intuition is that the latest support should remain valid in the near future, where the network has not changed much. FIFO ("First In First Out") defines the utility of a value by the time stamp when it was found or its validity was checked.

4.

Frequency. The expectation here is that the residue that has passed validity check frequently in the past will continue to do so in the future. FRQCYMIN initially sets the utility of a new support to one, then adds one to the score for every successful validity check for a residue and subtracts one for every unsuccessful check.

LEVELMIN, DOMMAX, FIFO and FRQCYMIN are called heuristic policies because they are based on the rationale for the corresponding factor. We also have the corresponding anti-heuristics which use the opposing rationale and whose anti-heuristic score is simply the negative of the heuristic score: LEVELMAX, DOMMAX, LIFO and FRQCYMAX.

In addition, we also consider a random replacement policy RANDOM as a baseline for comparison purpose in experiments.

Implementations of Different Policies

Intuitively, it makes sense that the residue store shouldn't be large -this is also borne out by our empirical study. Thus, we use an array for the residue store to speed up access. Other data structures may be more efficient with a larger residue store but have more overhead. Arrays also have better cache behavior [START_REF] Mitchell | A sat solver primer[END_REF]).

We use the following convention. Array indices range from 0 to maxR -1. Given an array D, D.size (initialized to 0) is the current number of residues. In the routine addResidue, D represents res[C, X, t[X]].

FIFO Policy

We use circular arrays to hold residues for static and fully dynamic FIFO. There is no need to record the utility score explicitly. Given an array D, D.head gives the index of the residue with the highest score; it is initially set to zero.

To find a support, both static and fully dynamic variants simply search from D.head until the end of the circular array (Algorithm 8). When a residue r is found to be valid in the circular array, fully dynamic FIFO sets D.head to point to r, since it has now become the most recent support. Moreover, the invalid residues before r are automatically placed at the end of the circular array as a side-effect. To add a value to the store, we simply store it before (in a modular way) D.head (Algorithm 7).

Algorithm 7: addResidue(C, X, t)

(FI F O) if D.size < maxR then D.size ← D.size + 1 1 if updateType = dynamic then 2 for p ← D.size -1 down to 1 do D[p] ← D[p -1] 3 D[0] ← t 4 else 5 D.head ← (D.head -1 + maxR) mod maxR 6 D[D.head] ← t 7
The residue store for dynamic FIFO also uses arrays. To find a support, we search array from beginning to the end. Since the utilities of the invalid residues before the first valid residue are not updated, when the first valid residue is found, we simply move it to the front and shift the invalid residues one position to the right (Algorithm 9). To add a new residue, we shift the first maxR -1 residues one position to the right and put the new one at the front (Algorithm 7). Array copying can be avoided but other alternative implementations may have higher overheads for the residue store sizes here. Here, only static and dynamic variants are described. The fully dynamic variant involves updating the utility of invalid residues checked so far as well as the utility of the first valid residue found. All these residues need to be re-ordered afterward. The cost of sorting these residues every time existValid is called is too expensive for our purpose. In contrast, the dynamic variant needs only to shift a single element in the array to its appropriate place to maintain the ordering of residues (lines 3-13 of Algorithm 11).

(dynamic FI F O) for p ← 0 to D.size -1 do 1 if isValid(C, D[p]) then 2 temp ← D[p] 3 for i ← p down to 1 do D[i] ← D[i -1] 4 D[0] ← temp
Algorithm 10: addResidue(C, X, t) with t/1000 denoting the probability of a tuple not allowed by a constraint. We used 21 instances for each class. The cpu time is the time (in seconds) it takes to solve 21 of them.

(static and dynamic LE V E L/DO M) if D.size < maxR then 1 D.size ← D.size + 1 2 else if D[D.size -1].score ≥ f (t) then return 3 p ← D.size -2 4 while p ≥ 0 ∧ D[p].score ≤ f (t) do 5 D[p + 1] ← D[p] 6 p ← p -1 7 D[p + 1] ← (t, f (t))
[p -1].score ≤ f (t) do 5 D[p] ← D[p -1] 6 p ← p -1 7 D[p] ← (t, f (t)) 8 if f (t) < D[p].
We also considered the following academic and real world problems: ehi-85-297-12, ehi-85-297-13, frb40-19-3, frb35-17-5, pigeons-10, pigeons-11, qa-5, qa-6, qk-20-0-5, qk-25-0-5, fapp01-0200-8, fapp01-0200-9, graph-10, graph-14, scen-11, scen-05 (called P1-P16).

For each problem, we collected the performance data using algorithms with different policies and residue store size from 1 to 10. Because it is obvious that an algorithm using larger number of residues requires fewer number of constraint checks, we made the cost of each constraint check as small as possible in order to test the performance in the worst case. Due to limited space, we can only present selected data, but all the analysis and observations apply to all the data unless mentioned otherwise.

HEURISTICS VERSUS ANTI-HEURISTICS

The experimental results on the effectiveness of heuristics and anti-heuristics are shown in Figure 1. The x-axis is an ordered sequence of pairs (problem, residueNumber) where problem is R1-R7 or P1-P16 and residueNumber is 1-10. We expect that the heuristics will always be better than the anti-heuristics. The results show that DOMMAX and LEVELMIN indeed performed better than DOMMIN and LEVELMAX on any data-point. The hypothesis is clearly supported by empirical evidence: the residue found when the domain size or the search level is smaller is more robust.

The result for FRQCY is surprising, however. This can be attributed to the fact that FRQCYMAX behaves in similar way to FIFO, which has been shown very effective and robust in our results. FRQCYMIN on the contrary favors new support less. This result indicates that the residues used most frequently might be less relevant in the future.

QUANTITATIVE EFFECT OF RESIDUE NUMBER ON PERFORMANCE

Having established that anti-heuristics are always worse, we will only consider the heuristics in the subsequent reports2 . We study in this section how the number of residues affects the broad performance of algorithms. Our data shows some clear patterns. First, the number of constraint checks decreases sharply and quickly converges to a stable number as the residue number increases. Second, the extra cost reflected mainly by validity checks increases at most linearly. As an illustration, we give performance of FIFO on all problems with varying residue number in Figure 2. The x-axis is the same as in the previous subsection. The two observations above imply that for big residue number, the new algorithms will not pay off. However, for small numbers (say 1 to 5), the extra cost can be compensated by the savings over constraint checks. The cpu time in Figure 2 reflects this interaction. The cpu time saving is not obvious due to the cheap cost of constraint checks. The third observation is that as the problems become more difficult, i.e., more constraint checks are needed, the savings (both in terms of constraint checks and cpu time) is more significant, as shown in the figure. The last observation is that the number of operations needed to maintain the proper ordering of the residue store is much lower than the number of the validity checks.

COMPARISON OF VARIOUS POLICIES AND UPDATING STRATEGIES

We compare performance of the static approach for FIFO, DOMMAX, LEVELMIN and FRQCYMIN on each problem class. Selected results on R7 and P15 are shown in Figure 3. Results on other problem classes show similar trends. For comparison purpose, we also tested a random replacement policy (RANDOM), described as follows. The residue store of RANDOM is an array of length maxR. The kth new residue found is inserted at the kth position in the array, where 1 ≤ k ≤ maxR. When the residue store is full (k > maxR), we pick an index within the range of the array randomly and overwrite the old residue with a new one. Searching for a support does not need to be undeterministic since the residue replacement already provides randomness for the policy. For simplicity, the search always starts from the lowest to the highest index in the array.

From the graphs, RANDOM is among the best policy for R7 while it is among the wost for P15. We notice that the relative performance of RANDOM also varies across different problem classes. Thus, the effect of RANDOM is dependent on many factors, which is not surprising. By contrast, the performance of other policies forms a clear and uniform ordering. Since its performance is unpredictable, we will no longer consider RANDOM in later analysis.

FIFO has the best results, followed by DOMMAX, LEVELMIN, and FRQCYMAX. DOM-MAX beats LEVELMIN possibly because it is finer-grained: LEVELMIN assigns the same utility to residues found at the same level regardless of individual differences in their respective domain size. While FRQCYMAX is closer to FIFO in the number of constraint checks, the cpu time is out of proportion with other policies. This is due partly to the cost of residue store ordering when the residue number gets larger, and partly to the large number of validity checks. From these graphs, we see that the ability of FRQCYMAX to retain supports is very close to that of FIFO, but the much larger number of validity checks implies that these supports are positioned at the very end of the array.

Next, we compare the performance of different updating strategies for FIFO in Figure 4. In cpu time, static FIFO is faster than dynamic FIFO when residue numbers are small, but becomes gradually slower as the number of residues increases. Nonetheless, the best results are obtained with static FIFO using low number of residues. The graphs for cpu time show that the best residue number lies in the region where the saving in constraint checks has just begun to be outweighed by the increase in the number of validity checks.

We do not give details for other policies but in general the difference between static and dynamic policy is small. For harder problems like P6 and P15, dynamic approach requires fewer number of constraint checks for DOMMAX and LEVELMIN while the converse is true for the anti-heuristics. However, as for cpu time, the dynamic approach is slower due to cheap constraint checks and higher cost in maintaining residue store dynamically.

Non-binary Problems

For non-binary problems, we used the Abscon solver. It implements MGAC embedding GAC3 rm with dom/wdeg variable ordering and lexicographical value ordering. Experiments were done in Linux on a cluster of 93 nodes, each with two Intel Xeon 3GHz and 2GB RAM. The results we present are for some classical hard non-binary instances from the Dimacs aim, Chessboard Coloration, Dubois, Schurr's Lemma, Dimacs Pret, Golomb's Ruler, All Interval Series, and Traveling Salesman Problem. Here, we have focused our attention to fully dynamic FIFO policy. 3 We only considered a limited number of residues (GAC3 rmk is GAC3 rm with k residues associated with each CN-value) as it appears to be the right approach. The results in the table further verify our observations on the impact of residue number over performance. Specifically, the number of constraint checks drops sharply and converges quickly. We observe that MGAC3 rm2 is a good compromise between saving constraint checks and improving cpu time.

Instances M GAC2001 M GAC3 M GAC3 rm M GAC3 rm2 M GAC3 rm3

Conclusion

We have generalized the existing work on single residue to multiple residues. We have a thorough investigation of the multiple residue approach including complexity analysis, the policies to manage the residue store, and an extensive empirical study of the effectiveness of the policies and the impact of store size on the performance. For heuristic policies, a key observation is that the number of constraint checks decreases quickly and converges to a stable number as the number of residues increases -indeed, in our experiments the harder the problem is, the larger the saving in the number of constraint checks. However, the extra cost reflected by validity checks also increases steadily as the number of residues increases. This suggests that the optimum number of residues should be small (say 1 to 5) since the total cost would be dominated by the cost of validity checks as more and more residues are in use.

  c 2008Lecoutre et al..

  ) in C. If C is tightness-bounded, then either c = O(1) and s = O(d), or c = O(d) and s = O(1) since c + s = d. It implies that the worst-case accumulated time complexity of seekSupport for a CN-value (C,X,a) is O(d + d) = O(d). The overall complexity of AC3 rm is then O(ed 2 ).

Algorithm 8 :

 8 existValid(C, D) (static and fully dynamic FI F O) for i ← 0 to D.size -1 do 1 p ← (D.head + i) mod D.size 2 if isValid(C, D[p]) then 3 if updateType = fullyDynamic then D.head ← p

  Like dynamic FIFO, we use standard arrays to hold residues. However, for LEVEL and DOM policies, given array D and index i, D[i] refers to the tuple (t, u) stored at the i th position, where t (D[i].residue) is the residue and u (D[i].score) is the utility. The residues are arranged in a descending order according to the utility, with ties broken by giving priority to the most recent support. addResidue and existValid are listed in Algorithm 10-11. LEVELMIN/MAX and DOMMIN/MAX are achieved by setting f as shown in Section 4.2.

8

  Algorithm 11: existValid(C, D) (static and dynamic LE V E L/DO M) for p ← 0 to D.size -1 do 1 if isValid(C, D[p].residue) then 2 if updateType= dynamic and f (t) = D[p].score then 3 if f (t) > D[p].score and p > 0 then 4 while p > 0 and D

  score and p < D.size -1 then 9 while p < D.size -1 and D[p + 1].score > f (t)

Figure 2 :Figure 4 :

 24 Figure 1: Comparisons of heuristics and anti-heuristics.

then return false 3 return true

4

Algorithm 4: seekSupport(C: Constraint, X: Variable, a: Value): Tuple

t ← ⊥ 1 while t = ⊤ do 2 if t ∈ rel(C) then return t 3 t ← setNextValid(C, X, a, t) 4 return ⊤ 5

Function addResidue(C: Constraint, X: Variable, t: Tuple) (not listed), assigns t to res[C, X

The exception is FRQCY, for which the anti-heuristic (FRQCYMAX) is better.

Unlike in the previous section, fully dynamic FIFO is better than both dynamic and static FIFO on these non-binary problems. We do not consider other policies for space reasons, as it is clear that FIFO is the best policy.

FRQCY Policy

We use arrays to hold residues in a descending order by their utility. It doesn't make sense to consider the static version since residues in the store would have a constant score of 1.

The function addResidue is shown in Algorithm 12. As the utility score for dynamic FRQ-CYMIN can only increase, it may be possible that all residues have a score greater than 1. Requiring the new support score to improve on the lowest score for fully dynamic FRQCYMIN would mean that the new support would never be recorded as its initial score is one. Rather, we replace the residue with the lowest score regardless of its value.

The function existValid is given in Algorithm 13. The principle is the following: when the residue at position p is found to be valid, its score is incremented, while, for fully dynamic FRQ-CYMIN, its score is decremented. When a residue is found to be valid, we have to re-order the array. It is performed by calling the routine re-order(D) (not described). Dynamic and fully dynamic FRQCYMAX are similar and have not been listed.

Experimental Results

We studied the performance of multiple residues with various store sizes and policies. The benchmarks are problems from the 2006 CSP solver competition 1 .

Binary Problems

For binary benchmarks, we used a solver written in C++ that implements MAC embedding AC3 r (i.e., AC3 with uni-directional residues). It includes dom/deg variable ordering and the lexicographical value ordering. The results were obtained on a DELL PowerEdge 1850 (two 3.6GHz Intel Xeon CPUs) in Linux.

We first considered some classes of binary random instances situated at the phase transition for satisfiability. These are (40, 8, 753, 100), (40, 11, 414, 200), (40, 16, 250, 350), (40, 25, 180, 500), (40, 40, 135, 650), (40, 80, 103, 800), (40, 180, 84, 900) (called R1-R7). A class is represented by (n, d, e, t) 1. See http://www.cril.univ-artois.fr/∼lecoutre/research/benchmarks/benchmarks.html