François Fillaux
email: francois.fillaux@upmc.fr

Alain Cousson

Neutron scattering studies of K 3 H(SO 4 ) 2 and K 3 D(SO 4 ) 2 : The particle-in-a-box model for the quantum phase transition

Keywords: numbers: 61.66.-f, 61.05.F-, 64.70.Tg Neutron diffraction, Hydrogen bonding, Quantum phase transition

In the crystal of K 3 H(SO 4 ) 2 , or K 3 D(SO 4 ) 2 , dimers (SO 4 )

are linked by strong centrosymmetric hydrogen or deuterium bonds whose O• • •O length is ≈ 2.50 Å. We address two open questions. (i) Are H or D sites split or not? (ii) Is there any structural counterpart to the phase transition observed for K 3 D(SO 4 ) 2 at T c ≈ 85.5 K, that does not exist for K 3 H(SO 4 ) 2 . Neutron diffraction by single-crystals at cryogenic or room temperature reveals no structural transition and no resolvable splitting of H or D sites. However, the width of the probability densities suggest unresolved splitting of the wavefunctions suggesting rigid entities H L1/2 -H R1/2 or D L1/2 -D R1/2 whose separation lengths are l H ≈ 0.16 Å or l D ≈ 0.25 Å. The vibrational eigenstates for the center of mass of H L1/2 -H R1/2 revealed by inelastic neutron scattering are amenable to a square-well and we suppose the same potential holds for D L1/2 -D R1/2 . In order to explain dielectric and calorimetric measurements of mixed crystals K 3 D (1-ρ) H ρ (SO 4 ) 2 (0 ≤ ρ ≤ 1), we replace the classical notion of order-disorder by the quantum notion of discernible (eg D L1/2 -D R1/2 ) or indiscernible (eg H L1/2 -H R1/2 ) components depending on the separation length of the split wavefunction. The discernible-indiscernible isostructural transition at finite temperatures is induced by a thermal pure quantum state or at 0 K by ρ.

This dependence on R OO accounts for temperature, pressure, and isotope effects. As R OO decreases, the potential operator for the stretching coordinate is expected to evolve from a double-well with a high barrier (compared to kT ) to a single-well for the shortest bonds. The latter is widely thought of as mimicking symmetric intermediates for proton transfer along the νOH coordinate in (bio)chemical reactions. This is an important field of investigations for computational chemistry and our purpose is to determine effective potentials based upon experiments for benchmark tests of theoretical methods.

In crystals, strong symmetric hydrogen bonds are encountered for small entities, such as chelates or dimers, well separated from each other. Bond lengths, symmetry and probability densities at proton sites are best characterized with single-crystal neutron diffraction (SXND), compared to x-ray (SXXD). Inelastic neutron scattering (INS) reveals OH transitions with spectacular contrast of intensity, compared to infrared or Raman. Probability densities in the momentum space can be determined with neutron Compton scattering (NCS).

The best characterized short symmetric hydrogen bond in a crystal is that of the chelate of the potassiumhydrogenmaleate, KH(OOC-CH=CH-COO), for which R OO = 2.427(1) Å at 20 K. 1 This is one of the shortest distance ever reported. INS spectra reveal numerous νOH transitions consistent with a funnel-shaped potential, such that the bare proton is localized at the center in the ground state, or largely delocalized over the oxygen orbitals in the excited states beyond ≈ 500 cm -1 . Similar potentials have been proposed for centrosymmetric dimers of potassiumhydrogenbistrifluoroacetate, KH(CF 3 COO) 2 and the cesium analogue, CsH(CF 3 COO) 2 , both with R OO = 2.436(4) Å at cryogenic temperatures. 2 In crystals, the existence of symmetric double-wells for longer bonds is controversial. For the crystal of potassiumdihydrogenphosphate (KDP, R OO ≈ 2.50 Å), Reiter et al. 3 claimed that NCS suggests a symmetric double-well in the paraelectric phase above T c ≈ 124 K, or a single minimum in the ferroelectric phase below T c . However, this is clearly excluded by the crystal symmetry and vibrational spectra. 4 Quasi-symmetric double wells have been reported for centrosymmetric dimers of the hydrogencarbonate family (MHCO 3 with M = K, Rb, Cs), 5,6 or benzoic acid, 7,8 for which R OO ≈ 2.60 Å, but these hydrogen bonds are not symmetric.

Apart from crystals, the chelate of the malonaldehyde (3-hydroxy-2-propenal) molecule in the gas phase (R OO ≈ 2.57 Å) is the best characterized symmetric double-well. 9,10 In the present work we report experimental and theoretical studies of the crystals of trispotassiumhydrogendisulfate, K 3 H(SO 4 ) 2 , and of the isostructural K 3 D(SO 4 ) 2 . [11][12][13] The unit cell of K 3 H(SO 4 ) 2 is comprised of 4 equivalent centrosymmetric dimers, SO 4

• • •H• • •SO 4
for which R OO ≈ 2.515(2) Å at 293 K or ≈ 2.487(2) Å at 100 K. These small dimers are well separate by K nuclei. The monoclinic group symmetry is A2/a and there is no structural phase-transition detectable by SXXD from 4 K to 300 K. Noda et al. 12 suggested that split proton sites at room temperature turn into single proton sites at low temperature, as R OO becomes shorter (geometric effect). INS reveals OH bending modes at 1250 or 1600 cm -1 and the νOH mode at ≈ 57 cm -1 . This is the lowest frequency ever reported for any hydrogen bond in a crystal. 14 This could be an evidence to suggest an asymmetric doublewell with a low barrier, but this is forbidden by the crystal symmetry. We report below SXND measurements showing that a double-well is unlikely at every temperature and we propose a new potential operator for INS.

Several members of the M 3 H(XO 4 ) 2 family (M = K, Rb, Cs; X = S, Se) demonstrate remarkable isotope effects which reveal intriguing differences between H and D bonds. In the specific case of K 3 D(SO 4 ) 2 , calorimetric and dielectric measurements reveal anomalies attributed to a phase transition at T c ≈ 85.5 K, whereas there is no transition for K 3 H(SO 4 ) 2 . 15 So far, T c is conceived of as an order-disorder transition from an antiferroelectric phase (A2) below T c to a paraelectric phase (A2/a) above T c , 15-17 A widespread explanation is that deuterons are delocalized in a symmetric double-well above T c and statistically localized below T c . 11,12,[17][18][19] Two models have been proposed to explain the isotope effect. (i)

The pseudospin Ising model with markedly different tunnel splitting for H or D. 15,17,19 (ii) Markedly different geometric effects for H or D bonds. 11,12,18,20 Both models are weakened by the lack of unquestionable evidence to suggest A2 symmetry consistent with asymmetric

O• • •D• • •O bonds below T c .
We have, therefore, carried out SXND measurements in order to observe whether or not there is a structural counterpart (possibly A2 ←→ A2/a) to the dielectric and calorimetric anomalies. This article is organized as follows. In Sec. II, we present SXND data for K 3 H(SO 4 ) 2 and K 3 D(SO 4 ) 2 at cryogenic and room temperatures. It transpires that the space group is A2/a at every temperature for both crystals and the best refined structures are obtained with H or D at the center of symmetry. Statistical distributions and antiferroelectricity below T c are, therefore, ruled out. In Sec. III, we show that SXND and INS are together consistent with unresolved split protons in a square-well and the same model is supposed to hold for deuterons. In Sec. IV we propose a new interpretation within the framework of quantum mechanics of the dielectric and calorimetric anomalies for KD (1-ρ) H ρ (SO 4 ) 2 , 0 ≤ ρ ≤ 1. We emphasize that these crystals are macroscopically quantum, to the least, below T c . In Sec.

V we show that the same model holds for the isostructural Rb 3 H(SO 4 ) 2 .

II. NEUTRON DIFFRACTION

Single-crystals were obtained by slow cooling of H

2 SO 4 (D 2 SO 4 );1.2 K 2 SO 4 in H 2 O(D 2 O).
Prismatic colorless specimens wrapped in aluminum were loaded in a closed-cycle-refrigerator whose temperature was controlled to ±1 K. Data were collected with the four-circle diffractometer 5C2 at the Orphée reactor (Laboratoire Léon Brillouin). 21 Every structural parameter computed with CRYSTALS 22 was allowed to vary independently. 23 Inspection of intensities for absent reflections confirms the monoclinic space group A2/a for both K 3 H(SO 4 ) 2 and K 3 D(SO 4 ) 2 at every temperature (see Table I). There is no visible symmetry breaking to suggest a statistical distribution of H or D. Best refinements were obtained with H or D at special positions (Tables II andIII) and every attempt to split these sites was less satisfactory. The best refined structures are comprised of Z = 4 dimer entities per unit cell. The slightly distorted SO 4 tetrahedra are consistent with spontaneous dielectric polarization. They are linked through centrosymmetric O

• • •H• • •O or O• • •D• • •O
bonds practically parallel to (a, b) and at ≈ ±30 • with respect to b (Fig. 1). The supposedly antiferroelectric phase of K 3 D(SO 4 ) 2 below T c is not confirmed and both the Ising model and geometric effects are rejected.

Unit cell parameters (Table I), positional parameters (Tables II andIII)), and thermal parameters for heavy nuclei (Tables IV andV the S-O bond lengths (R H SO ) are practically unchanged (Table VI). For K 3 D(SO 4 ) 2 , the unit cell is practically temperature independent, R ODO shortens very slightly from 2.532(1) Å to VII. The corresponding half-widths at half-height (HWHH), ∆ α (T ) (α = x, y, z), of the gaussian profiles [14]. The bands at ≈ 600 cm -1 were assigned to SO 4 entities.

P (α, T ) = exp [-α 2 /U α (T )] πU α (T ) , (1)

III. THE INFINITE SQUARE-WELL

Since the double-well proposed in Ref. [14] is rejected by SXND, we seek a singlewell potential for a proton oscillator (m = 1 amu) consistent with the νOH transition observed at (57 ± 3) cm -1 and with the probability density profile. A (quasi)harmonic model is inappropriate because the mean-square amplitude in the ground state, x 2 0 ≈ 0.29 Å2 , 24 should be one order of magnitude greater than U H

x (5). Consequently, we envisage another textbook potential that is the infinitely deep square-well such that V = 0 for [START_REF] Cohen-Tannoudji | Mécanique Quantique[END_REF] The eigenstates are:

-a H /2 ≤ x ≤ a H /2 and V = ∞ otherwise.
E n = n 2 h 2 8ma 2 H ψ n (x) = 2 a H sin nπ(x + a H /2) a H        ; n ≥ 1. (2)
The zero-point energy is

E H 1 = (E H 2 -E H 1 )/3 = (19 ± 1) cm -1 and a H = (0.47 ± 0.01) Å. Table IX shows that E H n -E H 1 = (n 2 -1)E H 1 , 2 ≤ n ≤ 6,
Ψ 2 (T ) = i exp(-(E i -E 1 )/kT )|ψ i | 2 i exp(-(E i -E 1 )/kT ) , (3)
can be confronted with (1). At 5 K, the leading term |ψ 1 | 2 gives a bell-shaped profile whose

HWHH is ∆ H 1x = a H /4 ≈ 0.
12 Å. This is significantly less than ∆ H x (5) ≈ 0.156 Å. We attribute the difference to a splitting of the wavefunction, what means that the energy levels refer to the center of mass of a rigid entity, say H L1/2 -H R1/2 , whose separation length is l H and whose probability density is the sum of two unresolved profiles |ψ 1 (x ± l H /2)| 2 /2. Our best estimate is l H = (0.16 ± 0.01) Å (see Fig. 4). Then, Ψ 2 (293) computed with the same parameters cannot be compared directly to a gaussian profile. The best we can conclude is that the HWHH corresponding to the plateau at half-height, namely (0.20 ± 0.05) Å, is in qualitative agreement with ∆ H

x (293). At this modest level of precision, there is no visible change of a H or l H consistent with the expansion of the O• • •O bond. There is no visible temperature induced geometric effect.

For K 3 D(SO 4 ) 2 , OD transitions are unknown. We tentatively suppose 14) is consistent with l D ≈ 0.25 Å (Fig. 5). (Note VII).

a D ≈ a H , E D n ≈ E H n /2, and |ψ D n | 2 ≈ |ψ H n | 2 . Then, ∆ D x (
that because of the mass effect, Ψ 2 (14) includes a minor contribution of

|ψ 2 (x ± l D /2)| 2 /2
to ≈ 5% of the total density). At 293 K, the HWHH of the calculated profile is roughly on the order of ∆ D x (293). Note that l Dl H ≈ 0.09 Å is rather close to R ODO -R OHO ≈ 0.07 Å at low temperatures. This suggests a geometric effect upon isotope substitution.

So far, the analysis of the density profiles is reasonably well established for K 3 H(SO 4 ) 2 at low temperature, since E H n and l H are determined by different measurements. For K 3 D(SO 4 ) 2 at low temperature, we ignore whether the assumption E D n ≈ E H n /2 is correct, so l D should be regarded with caution. In addition, it is impossible to conclude whether or not the model is pertinent at elevated temperatures. These drawbacks are wiped away in the next section where we show that dielectric and calorimetric data reported for mixed isotope crystals VII).

K 3 D (1-ρ) H ρ (SO 4 )

IV. ISOTOPE EFFECTS AND PHASE TRANSITIONS

Because centrosymmetric dimers have no electric dipole, there is no significant coupling terms between proton or deuteron vibrations which are, therefore, degenerate. This is confirmed by INS. The Ising model is, therefore, irrelevant.

Consider an ideal defect-free crystal of K 3 D(SO 4 ) 2 represented by a periodic lattice in 3-dimension of N indistinguishable rigid entities, D L1/2j -D R1/2j , centered at every node j.

N is on the order of Avogadro's constant. The centrosymmetric wavefunction for the center of mass x j is

ψ nj (x j , l D ) = 1 √ 2 [ψ n (x j + l D /2) -ψ n (x j -l D /2)], (4)
The lattice can be represented by N Bloch-wavefunctions indexed r (1 ≤ r ≤ N ),

Ψ n (x, l D , k nr ) = 1 N N j=1 ψ nj (x j , l D ) exp ik nr .L j , (5)
where x is a nonlocal observable independent of j, k nr is a wave vector and L j is a lattice (dǫ H /dT < 0 while dǫ D /dT > 0) and ǫ D reaches a maximum at T = T c . Above T c , probability densities, bond lengths and temperature effects for the dielectric constants are similar:

dǫ H /dT < 0; dǫ D /dT < 0. From these observations we tentatively infer that the S-O bond lengths, hence the dielectric constants, could correlate with the splitting of the wavefunction or, more specifically, with the nonlocal "indiscernibility degree" in the ground state proportional to the overlap defined on a 0 -1 scale as:

I 1 (l) = 2 a (a-l)/2 -(a-l)/2 sin π(x + a/2 + l/2) a sin π(x + a/2 -l/2) a dx = a -l a cos πl a + 1 π sin π a -l a            . (6)
Then, we tentatively set the discernibility-indiscernibility boundary (DIB) to I 1 (l di ) = 1/2, for l di = 0.183 Å. Consequently, I 1 (l H ) > 1/2 at every temperature (see Fig. 6 a) and the S-O bonds are virtually temperature independent. In contrast, I 1 (l D ) < 1/2 at 14 K and we assume that indiscernibility emerges at elevated temperatures from a thermal pure quantum state leading to coherent oscillations of the center of mass:

|Ψ D (t) = n≥1 α D n |ψ n exp(iE D n t/ ) |α D n | 2 = exp[-(E D n -E D 1 )/kT ] n≥1 exp[-(E D n -E D 1 )/kT ]          . (7)
The temperature law for the indiscernibility degree is then

I D (T ) = |α D 1 | 2 I 1 (l D ) + n>1 m>n β n,m |α D m | 2 , (8)
where β n,m is the indiscernibility degree of the superposition |n + e iϕ |m for which the mean position of the center of mass is: 25 and E D n ≈ E H n /2, our estimate l D = 0.25 Å and the thermal pure quantum state (7), are thus validated at the macroscopic level of the phase transition. We conclude that the square-well is largely isotope and temperature independent, to the least below T di .

x n,m (t) = X n,m cos E D m -E D n t + ϕ ; X n,m = 2 a a 0 x sin nπx a sin mπx a dx = a [(n -m)π] 2 [cos(n -m)π -1] - a [(n + m)π] 2 [cos(n + m)π -1]                  . ( 9
) If m = n+ 1, X n,m ≈ 0 and β n,m ≈ 0. Otherwise, X n,n+1 ≈ -2aπ -2
Our interpretation of the isotope effect is as follows. For K 3 H(SO 4 ) 2 , I 1 (l D ) > 1/2 at every temperature, the S-O bonds are practically temperature independent, and the monotonous decrease of ǫ H upon rising T is due to thermal fluctuations of the crystal lattice.

For K 3 D(SO 4 ) 2 , I 1 (l D ) < 1/2 for T < T di , the S-O bonds are longer and ǫ D < ǫ H . As T increases, I D (T ) increases, the S-O bonds shorten, ǫ D increases and reaches a maximum at T di = T c . Since the probability density computed at 85 K is similar to that at 293 K (see Fig. 5), we infer that the bond lengths are practically temperature independent above T di , so the monotonous decrease of ǫ D is analogous to that of ǫ H . So far, the cusp at T = T di for ǫ D reveals a discontinuity that is not accounted for by (8).

Consequently, we suppose that the phase transition is actually triggered by a discontinuity

IP = I 1 if l ≤ l di , IP = 1 if l ≥ l di .
of the indiscernibility parameter, say IP 1 , such that: (i) and β n,n+1 = 1 if T ≥ T di (see Fig. 7). This means that indiscernibility is quantified in the discernible domain whereas it is complete in the indiscernible range. This interpretation is corroborated by calorimetric measurements. For T < T di , we can label 1 or 2, respectively, the components of the split wavefunction (4) and

IP 1 = I 1 if l ≥ l di or IP 1 = 1 if l ≤ l di (see
ψ n (x j1 -l D /2) -ψ n (x j2 + l D /2) = ψ n (x j2 - l D /2) -ψ n (x j1 + l D /2). Otherwise, for T > T di , ψ n (x j1 -l D /2) -ψ n (x j2 + l D /2) ≡ ψ n (x j2 - l D /2) -ψ n (x j1 + l D /2).
Consequently, the total molar-entropy of the transition is R ln 2 and the specific heat jump is ∆C = 3R/2, as effectively measured. 15 The indiscernibility parameter is also relevant for K 3 D (1-ρ) H ρ (SO 4 ) 2 . Dielectric measurements realize Bloch-states at the center of the Brillouin-zone, via insignificant energy or momentum transfer, and variables of statistical-mechanical interest can be calculated from only one realization of the thermal pure quantum state. [START_REF] Sugiura | [END_REF] Protons and deuterons are not separate, so the indiscernibility parameter at 0 K is IP 1 (l ρ ), where

l ρ = ρl H + (1 -ρ)l D .
The DIB corresponds to a quantum phase transition at 0 K for ρ di = 0.74 ± 0.02. The

K 3 D (1-ρ) H ρ (SO 4 ) 2
, according to eqs ( 8) and (10). Dash: the dividing line between discernibility (below) and indiscernibility (above) domains.

temperature law analogous to (8),

IP ρ (T < T ρ di ) = |α D 1 | 2 IP 1 (l ρ ) + n>1 β|α D n | 2 IP ρ (T > T ρ di ) = 1    0 ≤ ρ < ρ di IP ρ (T ) = 1 ρ di < ρ ≤ 1          , (10)
accounts for the thermally induced phase transition at T ρ di for ρ < ρ di and for the lack of transition for ρ > ρ di . Here, we assume that E H n and E D n realized through energy transfer are separable and independent of ρ. Visual examination of Fig. 7 shows that for T < T ρ di , IP ρ is shifted upward by IP 1 (l ρ ) -IP 1 (l 0 ) (l 0 = l D ) and T ρ di can be written as:

T ρ di ≈ 1/2 -IP 1 (l ρ ) 1/2 -IP 1 (l 0 ) T 0 di ; 0 ≤ ρ < ρ di T ρ di = 0; ρ di < ρ ≤ 1      . ( 11
)
The quasi-linear variation of T ρ di for ρ < ρ di , mirroring that of IP 1 (l ρ ), is in reasonably good agreement with measurements (see Fig. 8). In contrast, the Ising model advocated by Moritomo et al. yields a curved function with a very steep slope as T c (ρ) -→ 0 (see Fig. 3 in Ref. [15]), so the extrapolated critical concentration (ρ c = 0.66 ± 0.04) 15 is smaller than ρ di . In practice, dielectric data suggest that it is technically difficult to measure T di <≈ 10 K with precision, so ρ c is logically a lower bound for ρ di . According to Fig. 7, the specific heat jump ∆C(ρ) for ρ < ρ di is proportional to

T ρ di 0 IP ρ (T )dT ≈ T ρ di [1/2 -IP 1 (l ρ )], (12)
so

∆C(ρ) ≈ 3R 2 × T ρ di [1/2 -IP 1 (l ρ )] T 0 di [1/2 -IP 1 (l 0 )] ; 0 ≤ ρ < ρ di ∆C(ρ) = 0; ρ di < ρ ≤ 1      . ( 13
)
This is in reasonably good agreement with observations (see Fig. 9). The shape of ∆R(ρ) mirrors that of IP 0 (T ) in Fig. 7. It is nearly linear for ρ <≈ 0.4 and the slope goes to zero as IP 1 (l ρ ) -→ IP 1 (l di ), so there is no discontinuity at the DIB.

V. Rb 3 H(SO 4 ) 2
The INS spectrum of the isomorphous crystal of Rb 3 H(SO 4 ) 2 (see Fig. 2 and Table VII) is also consistent with a square-well, apart from the transition observed at ≈ 44 cm -1 , instead of 57 cm -1 . However, the high background casts doubt on the actual frequency that should be reassessed with a more appropriate instrument. 14 In any case, an OH band in such a low frequency range is in conflict with the interpretation of NCS data proposed by Homouz et al. [START_REF] Homouz | Solid line: Temperature of the discernible-indiscernible transition for K 3 D (1-ρ) H ρ (SO 4 ) 2 , according to eq. (11)[END_REF] for a single-crystal of Rb 3 H(SO 4 ) 2 at 10 K. These authors observed that the mean = h/2a ≈ 1 Å-1 . [START_REF] Cohen-Tannoudji | Mécanique Quantique[END_REF] For the bending modes observed at ≈ 1250 cm -1 and ≈ 1600 cm -1 for both K 3 H(SO 4 ) 2 and Rb 3 H(SO 4 ) 2 , the harmonic approximation gives p 2 y 1/2 0 ≈ 4.39 Å-1 and p 2 z 1/2 0 ≈ 5.00 Å-1 .

Then, the projections along the crystal axes are:

p 2 b 1/2 = p 2 x 1/2 1 cos 30 • + p 2 z 1/2 0 sin 30 • ≈ 3.37 Å-1 p 2 a 1/2 = p 2 x 1/2 1 sin 30 • + p 2 z 1/2 0 cos 30 • ≈ 4.83 Å-1 p 2 c * 1/2 = p 2 y 1/2 0 ≈ 4.39 Å-1          . ( 14
)
These projections are effectively rather similar although the distribution is markedly anisotropic. Furthermore, the low νOH frequency imposes σ b < σ a , so we suspect σ a = 0.36 Å-1 ; leave open the question as to whether computational chemistry could confirm or not this counterintuitive well, as opposed to double-wells for longer bonds or funnel-shaped singlewells for shorter bonds.

σ b -p 2 a 1/2 1 = -0.23 Å-1 ; σ x -p 2 c * 1/2 1 = -0.04 Å-
Our explanation of the dielectric and calorimetric anomalies is based upon two postulates.

The first one is that these anomalies correspond to the crossing of a definite boundary for the discernibility-indiscernibility degree. To the best of our knowledge, this is unprecedented and we ignore whether or not this boundary is deeply rooted in the ground of quantum ref. [14]. The bands at ≈ 600 cm -1 were assigned to SO 4 entities. . . . . . . . . .

  Symmetric hydrogen bonds O• • •H• • •O are of fundamental importance to many disciplines across physics, chemistry, or biology, in the gaseous, liquid or solid state. Their properties depend on essentially the O• • •O bond length (say R OO ), such that a shorter R OO corresponds to a stronger bonding and a lower frequency for the OH stretching mode (νOH).

2 .

 2 520(2) Å (∆R ODO = -0.012(3) Å) and the R D SO 's increase significantly, what is in accordance with a change of the spontaneous dielectric polarization. The thermal parameters for H and D are markedly anisotropic. The diagonal factors U x (T ) and U z (T ) for coordinates x or z in (a, b) planes, parallel or perpendicular to O• • •O, respectively, are given in Table

FIG. 3 .

 3 FIG. 3. Eigenstates of the square-well along the OH stretching coordinate

FIG. 4 .

 4 FIG. 4. Solid lines: Computed probability densities along the OH stretching coordinate at 5 K (bottom) and 293 K (top) for a rigid H 1/2 -H 1/2 entity. The length is l H = 0.16 Å and the center of mass experiences a square-well whose width is 0.47 Å. Dashed line: probability density for each half proton. +: Gaussian profile whose variance is U Hx (T )/2 (see TableVII).

  2 are precisely represented by this model.

FIG. 5 .

 5 FIG. 5. Solid lines: Computed probability densities along the OD stretching coordinate at 14 K, 85 K, and 293 K for a rigid D L1/2 -D R1/2 entity. The length is l D = 0.25 Å and the center of mass experiences a box whose width is 0.47 Å. Dashed line: contributions of the n = 1 and n = 2 states to the probability density. +: Gaussian profile whose variance is U D x (T )/2 (TableVII).

  and the mean amplitude of oscillations, l = 2 √ 2aπ -2 , determines the time-averaged mean-separation of the splitwave: l D = l Dl = (0.11 ± 0.01) Å. Then, β n,n+1 = I 1 ( l D ) = 0.77 ± 0.03. The DIB I D (T di ) = 1/2 gives T di = (86±10) K, in accordance with T c . Our assumptions I 1 (l di ) = 1/2

FIG. 6 .

 6 FIG. 6. (a) The indiscernibility degree of the ground state, I 1 (l), as a function of the separation of the split-wave. l H = 0.16 Å, I 1 = 0.60 for K 3 H(SO 4 ) 2 . l D = 0.25 Å, I 1 = 0.27 for K 3 D(SO 4 ) 2 . l di = 0.183 Å, I 1di = 0.50 at the discernible-indiscernible boundary. (b) The transition parameter

Fig. 6 b

 6 ); (ii) IP D (T ) = I D (T ) and β n,n+1 = I 1 ( l D ), if T ≤ T di , or IP D (T ) = 1

FIG. 7 .

 7 FIG. 7. Temperature laws for the indiscernibility parameter driving the phase transition of

FIG. 8 .

 8 FIG. 8. Solid line: Temperature of the discernible-indiscernible transition for K3 D (1-ρ) H ρ (SO 4 ) 2 ,according to eq. (11). •: experimental data (dielectric measurements) estimated from digitized Fig.3in Ref.[15]. DIB: discernible-indiscernible boundary. l ρ = (1ρ)l D + ρl H is the effective separation length for the components of the wavefunction.

FIG. 9 . 2 x

 92 FIG.9. Solid: specific heat jump for K 3 D (1-ρ) H ρ (SO 4 ) 2 , according to eq. (13). Dot dash: the Ising model, as shown in Ref.[15]. •: experimental data estimated from digitized Fig.4in Ref.[15]. DIB: discernible-indiscernible boundary. l ρ = (1ρ)l D + ρl H is the effective separation length for the components of the wavefunction.

and σ b 2 b

 2 were swaped by Homouz et al. Then, taking into account this permutation, numerics (14) are in reasonably good agreement with NCS values: σ ap

  mechanics. The indiscernibility parameter can be thought of as a quantum counterpart to an order parameter diverging at the phase transition in classical physics. The second postulate is the correlation of the indiscernibility parameter and the dielectric constant, via the S-O bonds. The physical ground of this correlation is largely unknown. Our theory is scaled by three measured parameters, a, l H , l D . It accounts precisely for every observation, within the framework of the quantum theory of measurements. The crystal is a macroscopicscale quantum system for which the existence of a thermal pure quantum state is allowed by the adiabatic separation of H or D and heavy nuclei, as commonly encountered in similar O• • •O bonds. REFERENCES 1 F. Fillaux, N. Leygue, J. Tomkinson, A. Cousson, and W. Paulus, Chem. Phys. 244, 387 (1999).

  (T ) perpendicular to O• • •O and parallel to (a, b). U Hx (5) = 0.0351 U Hx (12) = 0.0326 U Hx (293) = 0.0932 U Dx (14) = 0.0937 U Dx (293) = 0.1117 U Hz (5) = 0.0162 U Hz (12) = 0.0160 U Hz (293) = 0.0290 U Dz (14) = 0.0141 U Dz (293) = 0.0254 TABLE VIII. Computed half-widths at half-height in Å units of the probability densities: ∆ H/Dx (T ) parallel to O• • •O; ∆ H/Dz (T ) perpendicular to O• • •O and parallel to (a, b); ∆ H/Dy (T ) parallel to c * . ∆ Hx (5) = 0.156 ∆ Hz (5) = 0.106 ∆ Hy (5) = 0.095 ∆ Hx (12) = 0.150 ∆ Hz (12) = 0.105 ∆ Hy (12) = 0.111 ∆ Hx (293) = 0.254 ∆ Hz (293) = 0.145 ∆ Hy (293) = 0.142 ∆ Dx (14) = 0.255 ∆ Dz (14) = 0.099 ∆ Dy (14) = 0.088 ∆ Dx (293) = 0.278 ∆ Dz (293) = 0.133 ∆ Dy (293) = 0.124

6 3 7 4 8 5 12 7

 67812 Eigenstates of the square-well along the OH stretching coordinate . . . . . . . . . . Solid lines: Computed probability densities along the OH stretching coordinate at 5 K (bottom) and 293 K (top) for a rigid H 1/2 -H 1/2 entity. The length is l H = 0.16 Å and the center of mass experiences a square-well whose width is 0.47 Å. Dashed line: probability density for each half proton. +: Gaussian profile whose variance is U H x (T )/2 (see Table VII). . . . . . . . . . . . . . . . . . . . . . . Solid lines: Computed probability densities along the OD stretching coordinate at 14 K, 85 K, and 293 K for a rigid D L1/2 -D R1/2 entity. The length is l D = 0.25 Å and the center of mass experiences a box whose width is 0.47 Å. Dashed line: contributions of the n = 1 and n = 2 states to the probability density. +: Gaussian profile whose variance is U D x (T )/2 (Table VII). . . . . . . . 9 6 (a) The indiscernibility degree of the ground state, I 1 (l), as a function of the separation of the split-wave. l H = 0.16 Å, I 1 = 0.60 for K 3 H(SO 4 ) 2 . l D = 0.25 Å, I 1 = 0.27 for K 3 D(SO 4 ) 2 . l di = 0.183 Å, I 1di = 0.50 at the discernible-indiscernible boundary. (b) The transition parameter IP = I 1 if l ≤ l di , IP = 1 if l ≥ l di . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Temperature laws for the indiscernibility parameter driving the phase transition of K 3 D (1-ρ) H ρ (SO 4 ) 2 , according to eqs (8) and (10). Dash: the dividing line between discernibility (below) and indiscernibility (above) domains. . . . . 13

  ) are similar for both isotopomers but temperature effects are different. For K 3 H(SO 4 ) 2 , the volume of the unit cell increases by ≈ 3% between 5

	5 K	14 K
	b	a
	293 K	293 K
	K 3 H(SO 4 ) 2	K 3 D(SO 4 ) 2
	FIG. 1. Projections of dimer entities onto (a, b) planes. Thermal ellipsoids represent 50% of the
	probability density for nuclei. Sulfur is yellow, oxygen is red, potassium is black, hydrogen or
	deuterium are grey.
	K and 293 K, R OHO increases from 2.463(1) Å to 2.496(1) Å (∆R OHO = 0.033(2) Å), whereas

  FIG. 2. INS spectra of K 3 H(SO 4 ) 2 and Rb 3 H(SO 4 ) 2 crystal powders at 20 K, after ref.

		0.3	57
	S(Q, ) (arbitrary units)	0.1 0.2	145	280	K 3 445	H(SO 4 600	) 2	680
		0	200		400		600	800	1000
									-1
				Energy Transfer (cm	)
		0.5	44
	S(Q, ) (arbitrary units)	0.1 0.2 0.3 0.4	145	260	Rb 3 440	H(SO 4 600	) 2	680
		0.0
		0	200		400		600	800	1000
									-1
				Energy Transfer (cm	)
	are given in Table VIII. Because the rather modest thermal parameters for heavy nuclei
	are virtually identical for both isotopomers, U α (T ) and ∆ α (T ) are largely representative
	of vibrational displacements. For K 3 H(SO 4 ) 2 , U H x (5) and U H x (293) suggest a markedly an-
	harmonic νOH mode at a low frequency and/or an unresolved splitting. For K 3 D(SO 4 ) 2 ,
	U D x (14) ≫ U H x (12) is at variance with the mass effect anticipated for usual single or double-
	well potentials.

  vector. Bloch-states are degenerate with respect to k nr , so the interpretation presented below is k nr independent. For K 3 H(SO 4 ) 2 , Bloch-states should be antisymmetrized if protons effectively behave as fermion, but this would be of no consequence to our interpretation.

It is worth emphasizing that a symmetric double-well giving rise to tunneling is excluded because there is no symmetry plane perpendicular to O• • •O compatible with a symmetric ground state 2 -1/2 [|L + |R ]. 26 Likewise, symmetric displacements of D L1/2j with respect to D R1/2j are forbidden. Only displacements of the center of mass of the rigid entity are allowed. Below T c , K 3 H(SO 4 ) 2 and K 3 D(SO 4 ) 2 demonstrate different probability densities for H or D, different S-O bond lengths, opposite temperature effects for their dielectric constants

  The space group symmetry of K 3 H(SO 4 ) 2 or K 3 D(SO 4 ) 2 determined with neutron diffraction is A2/a at cryogenic or room temperature and there is no resolvable site splitting for H or D. Statistical disorder and antiferroelectricity are ruled out. The anomalies of the dielectric constants or heat capacity do not correspond to any structural symmetry breaking. At low temperatures, the probability densities for H or D and the INS spectrum of K 3 H(SO 4 ) 2 are consistent with split-waves representing rigid entities H L1/2 -H R1/2 (l H ≈ 0.16 Å) or D L1/2 -D R1/2 (l D ≈ 0.25 Å) whose center of mass experiences a square-well. We

1 . The conflict of interpretation put forward by Homouz et al. is pointless and there is every reason to suppose that the square-well holds for Rb 3 H(SO 4 ) 2 . Our interpretation should be relevant for the dielectric anomaly of VI. CONCLUSION

TABLE I .

 I Neutron single crystal diffraction data and structure refinement for K 3 H(SO 4 ) 2 and K 3 D(SO 4 ) 2 . λ = 0.8305 Å. Space groups monoclinic A2/a with Z = 8. The criterion for used reflections was I > 3σ(I). The variance for the last digit is given in parentheses. Refinement on F.

		K 3 H(SO 4 ) 2			K 3 D(SO 4 ) 2
		5 K	12 K	293 K	14 K	293 K
	a ( Å)	9.704(1)	9.700(1)	9.806(1)	9.777(1)	9.778(1)
	b ( Å)	5.639(1)	5.634(1)	5.687(1)	5.674(1)	5.681(1)
	c ( Å)	14.526(1)	14.524(1)	14.702(1)	14.667(1)	14.701(1)
	β (˚)	102.95(1)	102.73(1)	102.93(1)	102.97(1)	103.05(1)
	V ( Å3 )	774.5(2)	774.3(2)	799.1(2)	792.9(2)	795.5(2)
	Measured reflections	2100	2385	2687	1822	794
	Independent reflections 1749	1373	1373	1253	432
	Used reflections	1524	1183	1781	999	404
	Coefficients	61	61	61	61	61
	R int	0.0002	0.0015	0.0002	0.0004	0.0006
	R-factor	0.027	0.045	0.024	0.047	0.018
	Weighted R-factor	0.023	0.051	0.014	0.033	0.020
	Goodness of fit	1.101	0.885	1.077	1.066	1.064
	Extinction coefficient	13.8(5)	14.8(10)	9.6(4)	5.0(3)	40(1)

TABLE II .

 II Atomic positions and isotropic temperature factors for K 3 H(SO 4 ) 2 at 5 K (first lines), 12 K (second lines) and 293 K (third lines).

	Atom	x/a	y/b	z/c	U(iso)( Å2 )
	K(1)	0.2500	0.75077(11)	0.0000	0.0041
		0.2500	0.7508(2)	0.0000	0.0048
		0.2500	0.75551(13)	0.0000	0.0198
	K(2)	0.34809(5)	0.73247(8)	0.30601(3)	0.0044
		0.34802(9)	0.73254(16)	0.30583(7)	0.0047
		0.34786(5)	0.73208(8)	0.30578(4)	0.0226
	S(1)	0.03807(6)	0.2726(1)	-0.11626(4)	0.0025
		0.03802(11)	0.2725(2)	-0.11641(8)	0.0024
		0.03807(5)	0.27302(8)	-0.11599(3)	0.0131
	O(1)	0.09896(3)	0.47808(5)	-0.15424(2)	0.0050
		0.09897(6)	0.4780(1)	-0.15423(5)	0.0055
		0.09819(3)	0.47694(5)	-0.15331(2)	0.0205
	O(2)	0.10868(3)	0.05211(5)	-0.13025(2)	0.0051
		0.10875(5)	0.0521(1)	-0.13021(4)	0.0053
		0.10805(4)	0.05494(5)	-0.13010(2)	0.0224
	O(3)	-0.11553(3)	0.25992(5)	-0.15300(2)	0.0050
		-0.11559(6)	0.2599(1)	-0.15298(4)	0.0052
		-0.11381(3)	0.26055(6)	-0.15244(2)	0.0205
	O(4)	0.06148(3)	0.31314(5)	-0.009127(19)	0.0052
		0.06141(6)	0.31321(11)	-0.00912(4)	0.0054
		0.06099(4)	0.31245(7)	-0.010210(19)	0.0241
	H(1)	0.0000	0.5000	0.0000	0.0213
		0.0000	0.5000	0.0000	0.0222
		0.0000	0.5000	0.0000	0.0477

TABLE III .

 III Atomic positions and isotropic temperature factors for K 3 D(SO 4 ) 2 at 14 K (first lines) 293 K (second lines).

	Atom	x/a	y/b	z/c	U(iso)( Å2 )
	K(1)	0.2500	0.7517(3)	0.0000	0.0042
		0.2500	0.7559(2)	0.0000	0.0217
	K(2)	0.34821(12)	0.7343(2)	0.30569(9)	0.0052
		0.34804(11)	0.73270(17)	0.30574(8)	0.0242
	S(1)	0.03846(15)	0.2704(3)	-0.1164(1)	0.0027
		0.03823(13)	0.2727(2)	-0.11604(8)	0.0142
	O(1)	0.09942(8)	0.47698(14)	-0.15382(5)	0.0048
		0.09827(6)	0.47659(11)	-0.15307(4)	0.0219
	O(2)	0.10852(9)	0.05043(15)	-0.13138(5)	0.0057
		0.10803(6)	0.05459(12)	-0.13039(5)	0.0238
	O(3)	-0.11496(8)	0.25884(15)	-0.15246(5)	0.0050
		-0.11348(7)	0.26004(11)	-0.15222(4)	0.0221
	O(4)	0.06275(8)	0.30902(15)	-0.00884(5)	0.0055
		0.06161(7)	0.31036(16)	-0.01022(4)	0.0261
	D(1)	0.0000	0.5000	0.0000	0.0396
		0.0000	0.5000	0.0000	0.0533

TABLE IV .

 IV Thermal parameters in Å2 units for K 3 H(SO 4 ) 2 at 5 K (first lines), 12 K (second lines) and 293 K (third lines). The thermal parameters account for the variation of the contribution of each atom to Bragg's peak intensities through the thermal factor T at depending on the reciprocal lattice parameters a * , b * , c * , and unit cell indexes in reciprocal space h, k, l, asT at = exp[-2π 2 (U at 11 h 2 a * 2 + U at 22 k 2 b * 2 + U at 33 l 2 c * 2 + 2Uat 12 ha * kb * + 2U at 23 kb * lc * + U at 31 lc * ha * )].

	Atom U 11	U 22	U 33	U 23	U 13	U 12
	K(1)	0.0043(2)	0.0048(2)	0.0031(2)	0.0000	0.00068(16)	0.0000
		0.0037(5)	0.0046(5)	0.0062(6)	0.0000	0.0010(4)	0.0000
		0.0184(2)	0.0242(2)	0.0165(2)	0.0000	0.00314(16)	0.0000
	K(2)	0.00446(16)	0.00402(16)	0.00458(15)	-0.00050(12)	0.00104(12)	-0.00024(12)
		0.0042(4)	0.0037(4)	0.0063(4)	-0.0008(3)	0.0010(3)	-0.0004(2)
		0.01899(15)	0.02015(16)	0.02811(19)	-0.00470(15)	0.00395(14)	-0.00035(14)
	S(1)	0.00272(19)	0.00217(18)	0.00252(18)	0.00001(14)	0.00057(14)	0.00014(15)
		0.0031(4)	0.0023(5)	0.0021(5)	-0.0004(3)	0.0009(3)	0.0001(3)
		0.01257(15)	0.01287(15)	0.01369(16)	0.00044(12)	0.00284(12)	0.00021(13)
	O(1)	0.00599(12)	0.00410(11)	0.00525(11)	0.00060(8)	0.00187(8)	-0.00108(8)
		0.0058(2)	0.0036(3)	0.0073(3)	0.00072(17)	0.00221(18)	-0.00122(15)
		0.02433(13)	0.0168(1)	0.02191(11)	0.00121(8)	0.00822(9)	-0.00460(8)
	O(2)	0.00578(11)	0.00396(11)	0.00583(11)	-0.00002(8)	0.00169(8)	0.00145(8)
		0.0054(2)	0.0033(3)	0.0074(3)	0.00026(17)	0.00179(19)	0.00153(15)
		0.02339(13)	0.01679(11)	0.02721(13)	0.00138(9)	0.00609(9)	0.00638(9)
	O(3)	0.00338(11)	0.00607(11)	0.0053(1)	-0.00023(8)	0.00057(7)	-0.00014(8)
		0.0034(2)	0.0056(3)	0.0063(3)	-0.00033(17)	0.00058(17)	-0.00018(14)
		0.01303(9)	0.02581(12)	0.02203(11)	-0.0023(1)	0.00271(7)	-0.00040(9)
	O(4)	0.00666(12)	0.00599(12)	0.0029(1)	-0.00011(7)	0.00090(8)	0.00021(8)
		0.0065(2)	0.0056(3)	0.0040(3)	-0.00022(17)	0.00078(17)	0.00022(15)
		0.02840(14)	0.03032(14)	0.0127(1)	0.00036(9)	0.00304(9)	0.0003(1)
	H(1)	0.0209	0.0304	0.0129	-0.0019	0.0041	-0.0082
		0.0204	0.0282	0.0178	-0.0017	0.0041	-0.0073
		0.0457	0.0765	0.0212	-0.0045	0.0081	-0.0282

TABLE V .

 V Thermal parameters in Å2 units for K 3 D(SO 4 ) 2 at 14 K (first lines) and 293 K (second lines). The thermal parameters account for the variation of the contribution of each atom to Bragg's peak intensities through the thermal factor T at depending on the reciprocal lattice parameters a * , b * , c * , and unit cell indexes in reciprocal space h, k, l, as T at = exp[-2π 2 (U at 11 h 2 a * 2 + U at 22 k 2 b * 2 + U at 33 l 2 c * 2 + 2U at 12 ha * kb * + 2U at 23 kb * lc * + U at 31 lc * ha * )].

	Atom	U 11	U 22	U 33	U 23	U 13	U 12
	K(1)	0.0037(6)	0.0056(7)	0.0032(6)	0.0000	0.0006(5)	0.0000
		0.0170(9)	0.028(1)	0.0196(8)	0.0000	0.0030(7)	0.0000
	K(2)	0.0043(4)	0.0041(5)	0.0071(4)	-0.0015(4)	0.0012(3)	-0.0003(4)
		0.0174(7)	0.0240(7)	0.0301(7)	-0.0057(5)	0.0033(5)	-0.0006(4)
	S(1)	0.0023(5)	0.0027(6)	0.0026(5)	0.0008(5)	-0.0002(4)	-0.0004(5)
		0.0109(7)	0.0155(7)	0.0161(7)	0.0019(5)	0.0026(5)	0.0004(4)
	O(1)	0.0057(3)	0.0030(3)	0.0060(3)	0.0006(2)	0.0020(2)	-0.0012(2)
		0.0227(4)	0.0200(5)	0.0243(5)	0.0010(3)	0.0079(3)	-0.0051(3)
	O(2)	0.0063(3)	0.0040(3)	0.0070(3)	0.0002(3)	0.0020(2)	0.0016(3)
		0.0216(4)	0.0205(5)	0.0293(4)	0.0011(3)	0.0059(3)	0.0072(3)
	O(3)	0.0033(3)	0.0058(3)	0.0059(3)	-0.0012(3)	0.0009(2)	-0.0004(3)
		0.0115(4)	0.0298(5)	0.0243(5)	-0.0028(3)	0.0029(3)	-0.0009(2)
	O(4)	0.0075(3)	0.0055(3)	0.0036(3)	-0.0013(3)	0.0016(2)	0.0018(3)
		0.0280(5)	0.0335(5)	0.0158(5)	0.0003(3)	0.0028(3)	0.0003(4)
	D(1)	0.035(1)	0.0728(16)	0.0112(6)	-0.0053(7)	0.0066(5)	-0.035(1)
		0.0464(7)	0.0907(11)	0.0222(6)	-0.0052(6)	0.0066(5)	-0.0370(8)

TABLE VI .

 VI Bond lengths in Å units of the SO 4 entities and their expansion ∆ with temperature.K 3 H(SO 4 ) 2 K 3 D(SO 4 ) 2

		5 K	293 K	∆R H SO	14 K	293 K	∆R D SO
	S1-O1	1.4636(6)	1.4623(5)	0.0013(11)	1.4754(17)	1.4585(14)	0.0169(31)
	S1-O2	1.4555(6)	1.4545(6)	0.0010(12)	1.4637(17)	1.4527(13)	0.0110(30)
	S1-O3	1.4680(6)	1.4674(5)	0.0006(11)	1.4758(16)	1.4602(15)	0.0156(31)
	S1-O4	1.5385(6)	1.5372(5)	0.0013(11)	1.5577(16)	1.5354(14)	0.0223(30)

TABLE VII .

 VII Diagonal thermal factors in Å2 units: U H/Dx (T ) parallel to O• • •O and U H/Dz

TABLE IX .

 IX Tentative assignment scheme of the OH stretching transitions observed with inelastic neutron scattering at 20 K, from Ref.[14], and comparison with the calculated transitions for a bare proton in an infinite square well whose width is a = 0.47 Å neutron scattering at 20 K, from Ref.[14], and comparison with the calculated transitions for a bare proton in an infinite square well whose width is a = 0.47 Å . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 LIST OF FIGURES 1 Projections of dimer entities onto (a, b) planes. Thermal ellipsoids represent 50% of the probability density for nuclei. Sulfur is yellow, oxygen is red, potassium is black, hydrogen or deuterium are grey. . . . . . . . . . . . . . . . . . . . . . 5 2 INS spectra of K 3 H(SO 4 ) 2 and Rb 3 H(SO 4 ) 2 crystal powders at 20 K, after

			1 -2	1 -3	1 -4	1 -5	1 -6
	K 3 H(SO 4 ) 2	Obs. (cm -1 )	57 ± 3	145 ± 5	280 ± 10	445 ± 15	680 ± 20
		Calc. (cm -1 )	57 ± 3	152 ± 8	285 ± 15	456 ± 24	665 ± 35
	Rb 3 H(SO 4 ) 2	Obs. (cm -1 )	44±?	145 ± 5	260 ± 10	440 ± 15	680 ± 20

inelastic

LIST OF TABLES

I

Neutron single crystal diffraction data and structure refinement for K 3 H(SO 4 ) 2 and K 3 D(SO 4 ) 2 . λ = 0.8305 Å. Space groups monoclinic A2/a with Z = 8.

The criterion for used reflections was I > 3σ(I).


Mireille El Gheche
email: el.gheche@univ-mlv.frpesquet@univ-mlv.fr

Jean-Christophe Pesquet

Joumana Farah
email: joumanafarah@usek.edu.lb

Mounir Kaaniche
email: kaaniche@telecom-paristech.frpesquet@telecom-paristech.fr

Béatrice Pesquet-Popescu

PROXIMAL SPLITTING METHODS FOR DEPTH ESTIMATION

Keywords: Stereo vision, disparity estimation, convex programming, variational methods, parallel proximal algorithm, proximity operator

Stereo matching is an active area of research in image processing. In a recent work, a convex programming approach was developed in order to generate a dense disparity field. In this paper, we address the same estimation problem and propose to solve it in a more general convex optimization framework based on proximal methods. More precisely, unlike previous works where the criterion must satisfy some restrictive conditions in order to be able to numerically solve the minimization problem, this work offers a great flexibility in the choice of the involved criterion. The method is validated in a stereo image coding framework, and the results demonstrate the good performance of the proposed parallel proximal algorithm.

INTRODUCTION

A stereo vision system generates two views, called left and right images, by recording two slightly different view angles of the same 3D scene. An important task in stereo vision is to find corresponding pixels which result from the projection of the same 3D point onto the two image planes. The displacement between the pixel coordinates of these points is called disparity. This information plays a crucial role in many application fields such as 3D reconstruction, intermediate view synthesis, obstacle detection, and stereo/multiview image compression [START_REF] Kaaniche | Dense disparity map representations for stereo image coding[END_REF]. The disparity estimation problem has been extensively studied in computer vision [START_REF] Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF]. Traditionally, disparity estimation algorithms are basically classified into two categories: local methods and global ones. Algorithms in the first category, where the disparity at each pixel depends only on intensity values within a local window, perform well in highly textured regions. However, they often produce noisy disparities in textureless regions and fail at occluded areas. These problems can be reduced by using global methods which aim at finding the disparity field that minimizes a global energy function over the entire image. For this purpose, several energy minimization algorithms have been proposed. The most common approaches are dynamic programming [START_REF] Veksler | Stereo correspondence by dynamic programming on a tree[END_REF], graph cuts [START_REF] Kolmogorov | Computing visual correspondence with occlusions using graph cuts[END_REF] and variational methods [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF]. While dynamic programming and graph cuts methods operate in a discrete manner, variational techniques work in a continuous space. Therefore, they possess the advantage of producing a disparity field with ideally infinite precision. Among these global approaches, it has been shown in [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF] that variational-based disparity estimation methods are among the most competitive techniques, for their preservation of depth discontinuities. In [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF], the dense disparity estimation problem has been formulated as a convex programming problem within a global variational approach. More precisely, a convex quadratic objective function has been employed and minimized subject to some appropriate convex constraints. The used optimization algorithm is able to handle a wide range of constraints. However, in order to satisfy the conditions of convergence of this algorithm, the cost function must be quadratic and strictly convex. In the present work, we reformulate the problem introduced in [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF] and we develop a parallel proximal algorithm which allows us to efficiently solve the problem over feasibility sets determined by multiple constraints modelling prior information. Our main contribution is to extend the approach in [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF] to a more general framework by relaxing the limitations on the choice of the similarity measure. The proposed algorithm offers a great flexibility in the choice of the convex cost function to be minimized, which, for example, may be based on the ℓ 1 -norm, an ℓ p -norm with p > 1 or the Kullback-Leibler divergence. The remainder of this paper is organized as follows. In Sec. 2, we present the stereo matching problem formulated as a constrained optimization problem. In Sec. 3, an efficient iterative solution is proposed based on a parallel proximal method. Experimental results, showing the accuracy of the resulting disparity maps as well as their applications to stereo image coding, are given in Sec. 4. Finally, some conclusions are drawn in Sec. 5.

PROBLEM STATEMENT

The process of stereo matching amounts to finding for each pixel in the left image I L a corresponding pixel in the right image I R . The matching problem is then equivalent to searching for the disparity field u which minimizes an error measure:

J(u) = (x,y)∈D φ(I L (x, y) -I R (x -u(x, y), y)) , (1)
where φ is assumed to be a proper lower-semicontinuous convex function from R to ] -∞, +∞] and D is the considered image domain. Similarity measures often used in literature are the Sum of Square Differences (SSD) and the Sum of Absolute Differences (SAD).

Cost function

Despite the convexity of φ, note that J is non-convex with respect to the field u. To circumvent this difficulty, we assume that an initial estimate ū of u is available and that the magnitude difference of the fields u and ū is small enough so that we can perform a Taylor expansion of the nonlinear term I R (xu, y) around an initial estimate ū as follows:

I R (x -u, y) ≃ I R (x -ū, y) -(u -ū) I x R (x -ū, y) , (2)
where

I x R (x -ū, y)
is the horizontal gradient of the disparity compensated right image. (For simplicity, we have dropped the dependence on (x, y) of u and u). Inserting linearization (2) in (1), we end up with the following convex cost function:

J(u) = (x,y)∈D\O φ(T (x, y) u -r(x, y)) (3)
where

T (x, y) = I x R (x -ū(x, y), y), (4) r(x, y) = I R (x -ū(x, y), y) + ū(x, y) T (x, y) -I L (x, y) (5)
and O denotes the occlusion area. Indeed, the occluded pixels, which correspond to points visible only in one image of the stereo pair need to be discarded [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF].

Optimizing this criterion is an ill-posed problem as the components of T may locally vanish. Thus, to solve the problem in a reliable manner, it is necessary to incorporate additional constraints. This can be performed by formulating the problem within a set theoretic framework [START_REF] Combettes | The foundations of set theoretic estimation[END_REF].

Convex constraints

The goal of set theoretic estimation is to obtain a feasible solution satisfying various constraints. Typically, each constraint is described by a closed convex set S i with i ∈ {1, ..., m}, in a Hilbert image space H. The intersection S of all the m sets, the feasibility set, constitutes the family of admissible solutions [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF]. In this context, the stereo matching problem can be formulated as a constrained optimization problem where the cost function given by ( 3) is minimized over the feasibility set. We have then to

Find u ∈ S = m i=1 S i such that J(u) = inf J(S). (6)
Constraint sets can generally be modelled as lower level sets:

∀i ∈ {1, ..., m},

S i = {u ∈ H | f i (u) ≤ δ i } , (7)
where

f i : H → R is a convex function and δ i ∈ R.
As mentioned before, the construction of convex constraints is derived from the properties of the field to be estimated. An example of possible prior knowledge is the range of the disparity values. We can restrict the variation of the values within the specific range [u min , u max ]. This can be expressed by the following constraint set S 1 :

S 1 = {u ∈ H | u min ≤ u ≤ u max } . ( 8
)
Most importantly, other constraints can be incorporated in order to enforce the smoothness of the disparity field in the homogeneous areas while preserving edges. First, we can use the Total Variation (TV) measure which attracted much attention in image recovery problems [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF]. Hence, a total variation based regularization constraint amounts to impose an upper bound τ > 0 on the TV:

S 2 = {u ∈ H | TV(u) ≤ τ } . (9)
Alternatively, we can adopt a wavelet domain approach to construct a smoothness constraint based on 2D separable wavelet coefficients. Inspired from the recent work in [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF], we consider the Besov space B 1 1,1 as being appropriate to model images containing discontinuities. The convex set associated with a semi-norm of this space is given by

S 3 = u ∈ H | j≥1,k∈Z 2 ,o∈{H,V } |c B j,k,o | ≤ κ (10)
where κ > 0 and the coefficients of u ∈ H in a wavelet basis or redundant frame are denoted by (c B j,k,o (u)) k∈Z 2 , o ∈ {H, V, D} being the orientation parameter and j ∈ N the resolution level.

PROPOSED ALGORITHM

Motivation

The objective of this section is to develop an iterative algorithm to solve the stereo matching problem as formulated by [START_REF] Combettes | The foundations of set theoretic estimation[END_REF]. The proximal method we propose allows us to find an optimal solution by using direct projections of the estimate onto the different constraint sets.

Optimization background

In Problem [START_REF] Combettes | The foundations of set theoretic estimation[END_REF], let us assume that H = R K and that each constraint set S i with i ∈ {1, . . . , m} can be expressed as L -1 i (C i ) where C i is a non-empty closed convex subset of R Ni and L i is a matrix in R Ni×K . Let us further assume that the projection onto each convex set C i takes a closed form (in the sense that it can be computed in a finite number of operations). The considered optimization problem then takes the following generic form: minimize Liu∈Ci,i∈{1,...,m}

J(u). ( 11
)
An algorithm allowing us to solve this problem is the algorithm described below. It basically consists of iterating computations of the proximity operator of J and projections onto convex sets. We recall that the proximity operator of the proper lowersemicontinuous convex function J is prox J : H → H : v → prox J v where prox J v is the unique minimizer of u → J(u) + 1 2 uv2 . This operator generalizes the notion of projection since prox ιC = P C where ι C is the indicator of a nonempty closed convex set C and P C denotes the projection onto C. For more details, the reader is referred to [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF].

PPXA+ algorithm

We propose to employ the PPXA+ (Parallel Proximal Algorithm) detailed in [START_REF] Pesquet | A parallel inertial proximal optimization method[END_REF], and which constitutes a generalization of some recent parallel proximal algorithms [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] Afonso | An augmented Lagrangian approach to the constrained optimization formulation of imaging inverse problems[END_REF]. As shown in Algorithm 1, it consists of computing parallel projections onto the sets (C i ) 1≤i≤m and the proximity operator of a scaled version of J, followed by a certain averaging procedure taking into account the matrices (L i ) 1≤i≤m . Suppose that the following assumptions hold.

• m i=1 L ⊤ i L i + I is an invertible matrix. • ∃ ǔ ∈ R K such that (∀i ∈ {1, . . . , m}) L i ǔ ∈ ri(C i ) and ǔ ∈ ri(dom J). 1 • (λ n ) n∈N is a sequence such that (∀n ∈ N) λ ≤ λ n+1 ≤ λ n < 2, where λ ∈]0, 2[.
Then, the sequence (u n ) n∈N generated by Algorithm 1 converges to a solution to Problem [START_REF] Combettes | The foundations of set theoretic estimation[END_REF] provided that such a solution exists.

SIMULATION RESULTS

Experiments were carried out on synthetic stereo images "corridor" and "teddy" taken from 2 .

[Initialization] (w 1 , ..., w

m ) ∈]0, +∞[ m ,γ > 0 (z i,0 ) 1≤i≤m+1 ∈ R N1 × . . . × R Nm × R K Q = ( m i=1 w i L ⊤ i L i + γI) -1 u 0 = Q( m i=1 w i L ⊤ i z i,0 + γz m+1,0 ) For n = 0, 1, . . . do For i = 1, . . . , m do p i,n = P Ci (z i,n ) end For p m+1,n = prox J γ (z m+1,n ) c n = Q( m i=1 w i L ⊤ i p i,n + γp m+1,n ) For i = 1, . . . , m do z i,n+1 = z i,n + λ n (L i (2c n -u n ) -p i,n ) end For z m+1,n+1 = z m+1,n + λ n (2c n -u n -p m+1,n ) u n+1 = u n + λ n (c n -u n ) end For Algorithm 1: Projective version of PPXA+.
In order to show the performance of the proposed method using an ℓ 1 error measure and the constraint sets S 1 and S 2 (here designated by "DDE-L1"), it has been compared with a block-based disparity estimation ("BDE") technique using the fixed size block matching method. For PPXA+ algorithm, we set w 1 = 100, w 2 = 10, γ = 200 and λ n ≡ 1.5. We have also tested the approach recently proposed by Miled et al. ("DDE-L2") which uses a cost function based on the ℓ 2norm and the same constraints. At this stage, it is worth pointing out that it was shown in [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF] that the quality of the results provided by the latter approach is comparable with state-ofthe-art methods, such as those based on graph cuts and dynamic programming [START_REF] Deng | A fast line segment based dense stereo algorithm using tree dynamic programming[END_REF]. Fig. 1 shows the generated disparity maps using the different mentioned methods. It can be seen that "DDE-L2" and "DDE-L1" produce smooth disparity maps. Furthermore, they are much less sensitive to noise than "BDE". Now, if we focus on the two dense disparity estimation methods, it can be observed that our approach "DDE-L1" better preserves the discontinuities around object boundaries. Numerical results are given by computing the PSNR values between the estimated disparity map and the ground truth. In order to emphasize the benefit from using the proposed method, we also evaluated its performance when used in a stereo image coding application. A straightforward approach consists of independently coding each image by using existing still image coders. However, most of the existing works rely on disparity compensation techniques [START_REF] Kaaniche | Dense disparity map representations for stereo image coding[END_REF]. The first step in these techniques consists of estimating the disparity map. Then, one image is considered as a reference image and the other is predicted in order to generate a prediction error referred to as a residual image. Finally, the disparity In order to alleviate this problem, we apply a quadtree coding of the estimated dense maps [START_REF] Kaaniche | Dense disparity map representations for stereo image coding[END_REF]. Fig. 2 illustrating the scalability in quality of the reconstruction procedure, shows that the proposed method outperforms "DDE-L2" by 0.1-0.3 dB.

Fig. 1 .

 1 Fig.1. Dense disparity maps for "corridor" and "teddy" stereo pair. From left to right: left image, "BDE" method, "DDE-L2" method, "DDE-L1" method .

Fig. 2 .

 2 Fig.2. PSNR (in dB) versus the bitrate (in bpp) after JPEG2000 encoding for the stereo pair "corridor".

dom J is the domain of J and the relative interior of a set C is denoted by riC.

http://cat.middlebury.edu/stereo/scenes2005/

CONCLUSIONS

We have presented a new parallel proximal algorithm for estimating disparity maps which has the advantage of allowing the use of various convex similarity measures and constraints. In particular, we have shown the benefit in coding applications which can be drawn from using an ℓ 1 cost function instead of the ℓ 2 one employed in previous DDE methods [5].


Ignacio Araya
email: iaraya@inf.utfsm.cl

Gilles Trombettoni
email: gilles.trombettoni@inria.fr

Bertrand Neveu
email: bertrand.neveu@enpc.fr

A Contractor Based on Convex Interval Taylor

Interval Taylor has been proposed in the sixties by the interval analysis community for relaxing continuous non-convex constraint systems. However, it generally produces a non-convex relaxation of the solution set. A simple way to build a convex polyhedral relaxation is to select a corner of the studied domain/box as expansion point of the interval Taylor form, instead of the usual midpoint. The idea has been proposed by Neumaier to produce a sharp range of a single function and by Lin and Stadtherr to handle n × n (square) systems of equations. This paper presents an interval Newton-like operator, called X-Newton, that iteratively calls this interval convexification based on an endpoint interval Taylor. This general-purpose contractor uses no preconditioning and can handle any system of equality and inequality constraints. It uses Hansen's variant to compute the interval Taylor form and uses two opposite corners of the domain for every constraint. The X-Newton operator can be rapidly encoded, and produces good speedups in constrained global optimization and constraint satisfaction. First experiments compare X-Newton with affine arithmetic.

Motivation

Interval B&B algorithms are used to solve continous constraint systems and to handle constrained global optimization problems in a reliable way, i.e., they provide an optimal solution and its cost with a bounded error or a proof of infeasibility. The functions taken into account may be non-convex and can include many (piecewise) differentiable operators like arithmetic operators (+, -, ., /), power, log, exp, sinus, etc.

Interval Newton is an operator often used by interval methods to contract/filter the search space [START_REF] Hansen | Global Optimization using Interval Analysis[END_REF]. The interval Newton operator uses an interval Taylor form to iteratively produce a linear system with interval coefficients. The main issue is that this system is not convex. Restricted to a single constraint, it forms a non-convex cone (a "butterfly"), as illustrated in Fig. 1-left. An n-dimensional constraint system is relaxed by an intersection of butterflies that is not convex either. (Examples can be found in [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF][START_REF] Kearfott | Rigorous Global Search: Continuous Problems[END_REF][START_REF] Moore | Introduction to Interval Analysis[END_REF].) Contracting optimally a box containing this non-convex relaxation has been proven to be NP-hard [START_REF] Kreinovich | Computational Complexity and Feasibility of Data Processing and Interval Computations[END_REF]. This explains why the interval analysis community has worked a lot on this problem for decades [START_REF] Hansen | Global Optimization using Interval Analysis[END_REF].

f' f' f' f' x x f(m([x]))+f'.(x-m([x])) f(m([x]))+f'.(x-m([x])) f(x)+f'.(x-x) f(x)+f'.(x-x) f(x) f(x)
Fig. 1. Relaxation of a function f over the real numbers by a function g : R → IR using an interval Taylor form (graph in gray). Left: Midpoint Taylor form, using a midpoint evaluation f (m([x])), the maximum derivative f ′ of f inside the interval [x] and the minimum derivative f ′ . Right: Extremal Taylor form, using an endpoint evaluation f (x), f ′ and f ′ .

Only a few polynomial time solvable subclasses have been studied. The most interesting one has been first described by Oettli and Prager in the sixties [START_REF] Oettli | On the Solution Set of a Linear System with Inaccurate Coefficients[END_REF] and occurs when the variables are all non-negative or non-positive. Unfortunately, when the Taylor expansion point is chosen strictly inside the domain (the midpoint typically), the studied box must be previously split into 2 n subproblems/quadrants before falling in this interesting subclass [START_REF] Aberth | The Solution of Linear Interval Equations by a Linear Programming Method[END_REF][START_REF] Beaumont | Algorithmique pour les intervalles[END_REF][START_REF] Chabert | Techniques d'intervalles pour la résolution de systèmes d'intervalles[END_REF]. Hansen and Bliek independently proposed a sophisticated and beautiful algorithm for avoiding explicitly handling the 2 n quadrants [START_REF] Hansen | Bounding the Solution of Interval Linear Equations[END_REF][START_REF] Bliek | Computer Methods for Design Automation[END_REF]. However, the method is restricted to n × n (square) systems of equations (no inequalities). Also, the method requires the system be first preconditioned (i.e., the interval Jacobian matrix must be multiplied by the inverse matrix of the domain midpoint). The preconditioning has a cubic time complexity, implies an overestimate of the relaxation and requires non-singularity conditions often met only on small domains, at the bottom of the search tree.

In 2004, Lin & Stadtherr [START_REF] Lin | LP Strategy for the Interval-Newton Method in Deterministic Global Optimization[END_REF] proposed to select a corner of the studied box, instead of the usual midpoint. Graphically, it produces a convex cone, as shown in Fig. 1-right. The main drawback of this extremal interval Taylor form is that it leads to a larger system relaxation surface. The main virtue is that the solution set belongs to a unique quadrant and is convex. It is a polytope that can be (box) hulled in polynomial-time by a linear programming (LP) solver: two calls to an LP solver compute the minimum and maximum values in this polytope for each of the n variables (see Section 4). Upon this extremal interval Taylor, they have built an interval Newton restricted to square n × n systems of equations for which they had proposed in a previous work a specific preconditioning. They have presented a corner selection heuristic optimizing their preconditioning. The selected corner is common to all the constraints.

The idea of selecting a corner as Taylor expansion point is mentioned, in dimension 1, by A. Neumaier (see page 60 and Fig. 2.1 in [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF]) for computing a range enclosure (see Def. 1) of a univariate function. Neumaier calls this the linear boundary value form. The idea has been exploited by Messine and Laganouelle for lower bounding the objective function in a Branch & Bound algorithm for unconstrained global optimization [START_REF] Messine | Enclosure Methods for Multivariate Differentiable Functions and Application to Global Optimization[END_REF].

McAllester et al. also mention this idea in [START_REF] Mcallester | Three Cuts for Accelerated Interval Propagation[END_REF] (end of page 2) for finding cuts of the box in constraint systems. At page 211 of Neumaier's book [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF], the step (4) of the presented pseudo-code also uses an endpoint interval Taylor form for contracting a system of equations. 1

Contributions

We present in this paper a new contractor, called X-Newton (for eXtremal interval Newton), that iteratively achieves an interval Taylor form on a corner of the studied domain. X-Newton does not require the system be preconditioned and can thus reduce the domains higher in the search tree. It can treat well-constrained systems as well as under-constrained ones (with fewer equations than variables and with inequalities), as encountered in constrained global optimization. The only limit is that the domain must be bounded, although the considered intervals, i.e., the initial search space, can be very large.

This paper experimentally shows that such a contractor is crucial in constrained global optimization and is also useful in continuous constraint satisfaction where it makes the whole solving strategy more robust.

After the background introduced in the next section, we show in Section 3 that the choice of the best expansion corner for any constraint is an NP-hard problem and propose a simple selection policy choosing two opposite corners of the box. Tighter interval partial derivatives are also produced by Hansen's recursive variant of interval Taylor. Section 4 details the extremal interval Newton operator that iteratively computes a convex interval Taylor form. Section 5 highlights the benefits of X-Newton in satisfaction and constrained global optimization problems.

This work provides an alternative to the two existing reliable (interval) convexification methods used in global optimization. The Quad [START_REF] Lebbah | Efficient and safe global constraints for handling numerical constraint systems[END_REF][START_REF] Lebbah | An Efficient and Safe Framework for Solving Optimization Problems[END_REF] method is an interval reformulation-linearization technique that produces a convex polyhedral approximation of the quadratic terms in the constraints. Affine arithmetic produces a polytope by replacing in the constraint expressions every basic operator by specific affine forms [START_REF] De Figueiredo | Affine Arithmetic: Concepts and Applications[END_REF][START_REF] Vu | Enhancing Numerical Constraint Propagation using Multiple Inclusion Representations[END_REF][START_REF] Baharev | Computation of an Extractive Distillition Column with Affine Arithmetic[END_REF]. It has been recently implemented in an efficient interval B&B [START_REF] Ninin | A Reliable Affine Relaxation Method for Global Optimization[END_REF]. Experiments provide a first comparison between this affine arithmetic and the corner-based interval Taylor.

Background

Intervals allow reliable computations on computers by managing floating-point bounds and outward rounding.

Intervals

An interval [x i ] = [x i , x i ] defines the set of reals x i s.t. x i ≤ x i ≤ x i , where x i and x i are floating-point numbers. IR denotes the set of all intervals. The size or

width of [x i ] is w([x i ]) = x i -x i . A box [x] is the Cartesian product of intervals [x 1 ] × ... × [x i ] × ... × [x n ]. Its width is defined by max i w([x i ]). m([x]) denotes the middle of [x]. The hull of a subset S of R n is the smallest n-dimensional box enclosing S.
Interval arithmetic [START_REF] Moore | Interval Analysis[END_REF] has been defined to extend to IR elementary functions over R. For instance, the interval sum is defined by

[x 1 ]+[x 2 ] = [x 1 +x 2 , x 1 +x 2 ].
When a function f is a composition of elementary functions, an extension of f to intervals must be defined to ensure a conservative image computation.

Definition 1 (Extension of a function to IR; inclusion function; range enclosure) Consider a function f : R n → R.

[f ] : IR n → IR is said to be an extension of f to intervals iff:

∀[x] ∈ IR n [f ]([x]) ⊇ {f (x), x ∈ [x]} ∀x ∈ R n f (x) = [f ](x)
The natural extension [f ] N of a real function f corresponds to the mapping of f to intervals using interval arithmetic. The outer and inner interval linearizations proposed in this paper are related to the first-order interval Taylor extension [START_REF] Moore | Interval Analysis[END_REF], defined as follows:

[f ] T ([x]) = f ( ẋ) + i [a i ] . ([x i ] -ẋi )
where ẋ denotes any point in [x], e.g., m([x]), and

[a i ] denotes ∂f ∂xi N ([x]).
Equivalently, we have:

∀x ∈ [x], [f ] T ([x]) ≤ f (x) ≤ [f ] T ([x]). Example. Consider f (x 1 , x 2 ) = 3x 2 1 + x 2 2 + x 1 x 2 in the box [x] = [-1, 3] × [-1, 5]. The natural evaluation provides: [f ] N ([x 1 ], [x 2 ]) = 3[-1, 3] 2 + [-1, 5] 2 + [-1, 3][-1, 5] = [0, 27] + [0, 25] + [-5, 15] = [-5, 67]. The partial derivatives are: ∂f ∂x1 (x 1 , x 2 ) = 6x 1 +x 2 , [ ∂f ∂x1 ] N ([-1, 3], [-1, 5]) = [-7, 23], ∂f ∂x2 (x 1 , x 2 ) = x 1 +2x 2 , [ ∂f ∂x2 ] N ([x 1 ], [x 2 ]) = [-3, 13]. The interval Taylor evaluation with ẋ = m([x]) = (1, 2) yields: [f ] T ([x 1 ], [x 2 ]) = 9 + [-7, 23][-2, 2] + [-3, 13][-3, 3] = [-76, 94].

A simple convexification based on interval Taylor

Consider a function f : R n → R defined on a domain [x], and the inequality constraint f (x) ≤ 0. For any variable

x i ∈ x, let us denote [a i ] the interval partial derivative ∂f ∂xi N ([x]
). The first idea is to lower tighten f (x) with one of the following interval linear forms that hold for all x in [x].

f (x) + a 1 y l 1 + ... + a n y l n ≤ f (x) (1)
f (x) + a 1 y r 1 + ... + a n y r n ≤ f (x) (2)
where: [START_REF] Araya | Exploiting Monotonicity in Interval Constraint Propagation[END_REF]. When applied to a set of inequality and equality2 constraints, we obtain a polytope enclosing the solution set.

y l i = x i -x i and y r i = x i -x i . A corner of the box is chosen: x in form (1) or x in form
The correctness of relation ( 1) -see for instance [START_REF] Trombettoni | Inner Regions and Interval Linearizations for Global Optimization[END_REF][START_REF] Lin | LP Strategy for the Interval-Newton Method in Deterministic Global Optimization[END_REF] -lies on the simple fact that any variable

y l i is non-negative since its domain is [0, d i ], with d i = w([y l i ]) = w([x i ]) = x i -x i . Therefore, minimizing each term [a i ] y l i for any point y l i ∈ [0, d i ] is obtained with a i . Symmetrically, relation (2) is correct since y r i ∈ [-d i , 0] ≤ 0,
and the minimal value of a term is obtained with a i . Note that, even though the polytope computation is safe, the floating-point round-off errors made by the LP solver could render the hull of the polytope unsafe. A cheap post-processing proposed in [START_REF] Neumaier | Safe Bounds in Linear and Mixed-Integer Programming[END_REF], using interval arithmetic, is added to guarantee that no solution is lost by the Simplex algorithm.

3 Extremal interval Taylor form

Corner selection for a tight convexification

Relations (1) and ( 2) consider two specific corners of the box [x]. We can remark that every other corner of [x] is also suitable. In other terms, for every variable x i , we can indifferently select one of both bounds of [x i ] and combine them in a combinatorial way: either x i in a term a i (x i -x i ), like in relation [START_REF] Aberth | The Solution of Linear Interval Equations by a Linear Programming Method[END_REF], or x i in a term a i (x i -x i ), like in relation [START_REF] Araya | Exploiting Monotonicity in Interval Constraint Propagation[END_REF].

A natural question then arises: Which corner x c of [x] among the 2 n -set X c ones produces the tightest convexification? If we consider an inequality f (x) ≤ 0, we want to compute a hyperplane f l (x) that approximates the function, i.e., for all x in [x] we want:

f l (x) ≤ f (x) ≤ 0.
Following the standard policy of linearization methods, for every inequality constraint, we want to select a corner x c whose corresponding hyperplane is the closest to the non-convex solution set, i.e., adds the smallest volume. This is exactly what represents Expression (3) that maximizes the Taylor form for all the points x = {x 1 , ..., x n } ∈ [x] and adds their different contributions: one wants to select a corner x c from the set of corners X c such that:

max x c ∈X c x1 x1=x1 ... xn xn=xn (f (x c ) + i z i ) dx n ... dx 1 (3)
where:

z i = a i (x i -x i ) iff x c i = x i , and z i = a i (x i -x i ) iff x c i = x i . Since:
f (x c ) is independent from the x i values, -any point z i depends on x i but does not depend on x j (with

j = i), - xi xi=xi a i (x i -x i )dx i = a i di yi=0 y i dy i = a i 0.5 d 2 i , - xi xi=xi a i (x i -x i )dx i = a i 0 -di y i dy i = -0.5 a i d 2 i , Expression (3)
is equal to:

max x c ∈X c i d i f (x c ) + i d i i 0.5 a c i d i
where

d i = w([x i ]
) and a c i = a i or a c i = -a i . We simplify by the positive factor i d i and obtain:

max x c ∈X c f (x c ) + 0.5 i a c i d i (4)
Unfortunately, we have proven that this maximization problem (4) is NPhard.

Proposition 1 (Corner selection is NP-hard) Consider a polynomial3 f : R n → R, with rational coefficients, and defined on a domain [x] = [0, 1] n . Let X c be the 2 n -set of corners, i.e., in which every component is a bound 0 or 1. Then,

max x c ∈X c -(f (x c ) + 0.5 i a c i d i ) (or min x c ∈X c f (x c ) + 0.5 i a c i d i ) is an NP-hard problem.
The extended paper [START_REF] Araya | A Contractor Based on Convex Interval Taylor[END_REF] shows straightforward proofs that maximizing the first term of Expression 4 (f (x c )) is NP-hard and maximizing the second term 0.5 i a c i d i is easy, by selecting the maximum value among a i and -a i in every term. However, proving Proposition 1 is not trivial (see [START_REF] Araya | A Contractor Based on Convex Interval Taylor[END_REF]) and has been achieved with a polynomial reduction from a subclass of 3SAT, called BALANCED-3SAT. 4Even more annoying is that experiments presented in Section 5 suggest that the criterion (4) is not relevant in practice. Indeed, even if the best corner was chosen (by an oracle), the gain in box contraction brought by this strategy w.r.t. a random choice of corner would be not significant. This renders pointless the search for an efficient and fast corner selection heuristic.

This study suggests that this criterion is not relevant and leads to explore another criterion. We should notice that when a hyperplane built by endpoint interval Taylor removes some inconsistent parts from the box, the inconsistent subspace more often includes the selected corner x c because the approximation at this point is exact. However, the corresponding criterion includes terms mixing variables coming from all the dimensions simultaneously, and makes difficult the design of an efficient corner selection heuristic.

This qualitative analysis nevertheless provides us rationale to adopt the following policy.

Using two opposite corners

To obtain a better contraction, it is also possible to produce several, i.e., c, linear expressions lower tightening a given constraint f (x) ≤ 0. Applied to the whole system with m inequalities, the obtained polytope corresponds to the intersection of these c m half-spaces. Experiments (see Section 5.2) suggest that generating two hyperplanes (using two corners) yields a good ratio between contraction (gain) and number of hyperplanes (cost). Also, choosing opposite corners tends to minimize the redundancy between hyperplanes since the hyperplanes remove from the box preferably the search subspaces around the selected corners.

Note that, for managing several corners simultaneously, an expanded form must be adopted to put the whole linear system in the form Ax-b before running the Simplex algorithm. For instance, if we want to lower tighten a function f (x) by expressions (1) and (2) simultaneously, we must rewrite:

1. f (x) + i a i (x i -x i ) = f (x) + i a i x i -a i x i = i a i x i + f (x) -i a i x i 2. f (x) + i a i (x i -x i ) = f (x) + i a i x i -a i x i = i a i x i + f (x) -i a i x i
Also note that, to remain safe, the computation of constant terms a i x i (resp. a i x i ) must be achieved with degenerate intervals:

[a i , a i ] [x i , x i ] (resp. [a i , a i ] [x i , x i ]).

Preliminary interval linearization

Recall that the linear forms (1) and ( 2) proposed by Neumaier and Lin & Stadtherr use the bounds of the interval gradient, given by ∀i ∈ {1, ..., n},

[a i ] = ∂f ∂xi N ([x]
). Eldon Hansen proposed in 1968 a variant in which the Taylor form is achieved recursively, one variable after the other [START_REF] Hansen | On Solving Systems of Equations Using Interval Arithmetic[END_REF][START_REF] Hansen | Global Optimization using Interval Analysis[END_REF]. The variant amounts in producing the following tighter interval coefficients:

∀i ∈ {1,..., n}, [a i ] = ∂f ∂x i N ([x 1 ] × ... × [x i ] × ẋi+1 × ... × ẋn )
where ẋj ∈ [x j ], e.g., ẋj = m([x j ]). By following Hansen's recursive principle, we can produce Hansen's variant of the form (1), for instance, in which the scalar coefficients a i are:

∀i ∈ {1, ..., n}, a i = ∂f ∂x i N ([x 1 ]×...×[x i ]×x i+1 ×...×x n ).
We end up with an X-Taylor algorithm (X-Taylor stands for eXtremal interval Taylor) producing 2 linear expressions lower tightening a given function f : R n → R on a given domain [x]. The first corner is randomly selected, the second one is opposite to the first one.

eXtremal interval Newton

We first describe in Section 4.1 an algorithm for computing the (box) hull of the polytope produced by X-Taylor. We then detail in Section 4.2 how this X-NewIter procedure is iteratively called in the X-Newton algorithm until a quasi-fixpoint is reached in terms of contraction.

X-Newton iteration

Algorithm 1 describes a well-known algorithm used in several solvers (see for instance [START_REF] Lebbah | Efficient and safe global constraints for handling numerical constraint systems[END_REF][START_REF] Baharev | Computation of an Extractive Distillition Column with Affine Arithmetic[END_REF]). A specificity here is the use of a corner-based interval Taylor form (X-Taylor) for computing the polytope. All the constraints appear as inequality constraints f j (x) ≤ 0 in the vector/set f = (f 1 , ..., f j , ..., f m ). x = (x 1 , ..., x i , ..., x n ) denotes the set of variables with domains [x].

The first loop on the constraints builds the polytope while the second loop on the variables contracts the domains, without loss of solution, by calling a Simplex algorithm twice per variable. When embedded in an interval B&B for constrained global optimization, X-NewIter is modified to also compute a lower bound of the objective in the current box: an additional call to the Simplex algorithm minimizes an X-Taylor relaxation of the objective on the same polytope.

Heuristics mentioned in [START_REF] Baharev | Computation of an Extractive Distillition Column with Affine Arithmetic[END_REF] indicate in which order the variables can be handled, thus avoiding in practice to call 2n times the Simplex algorithm.

X-Newton

The procedure X-NewIter allows one to build the X-Newton operator (see Algorithm 2). Consider first the basic variant in which CP-contractor = ⊥.

Algorithm 2 X-Newton (f , x, [x], ratio fp, CP-contractor): [x] repeat [x]save ← [x] [x] ← X-NewIter (f , x, [x]) if CP-contractor = ⊥ and gain([x],[x]save) > 0 then [x] ← CP-contractor(f ,x,[x]) end if until empty([x]) or gain([x],[x]save) < ratio fp) return [x]
X-NewIter is iteratively run until a quasi fixed-point is reached in terms of contraction. More precisely, ratio fp is a user-defined percentage of the interval size and:

gain([x ′ ], [x]) := max i w([x i ]) -w([x ′ i ]) w([x i ])
.

We also permit the use of a contraction algorithm, typically issued from constraint programming, inside the main loop. For instance, if the user specifies CP-contractor=Mohc and if X-NewIter reduces the domain, then the Mohc constraint propagation algorithm [START_REF] Araya | Exploiting Monotonicity in Interval Constraint Propagation[END_REF] can further contract the box, before waiting for the next choice point. The guard gain([x], [x] save ) > 0 guarantees that CP-contractor will not be called twice if X-NewIter does not contract the box.

Remark

Compared to a standard interval Newton, a drawback of X-Newton is the loss of quadratic convergence when the current box belongs to a convergence basin. It is however possible to switch from an endpoint Taylor form to a midpoint one and thus be able to obtain quadratic convergence, as detailed in [START_REF] Araya | A Contractor Based on Convex Interval Taylor[END_REF].

Also note that X-Newton does not require the system be preconditioned so that this contractor can cut branches early during the tree search (see Section 5.2). In this sense, it is closer to a reliable convexification method like Quad [START_REF] Lebbah | Efficient and safe global constraints for handling numerical constraint systems[END_REF][START_REF] Lebbah | An Efficient and Safe Framework for Solving Optimization Problems[END_REF] or affine arithmetic [START_REF] Ninin | A Reliable Affine Relaxation Method for Global Optimization[END_REF].

Experiments

We have applied X-Newton to constrained global optimization and to constraint satisfaction problems.

Experiments in constrained global optimization

We have selected a sample of global optimization systems among those tested by Ninin et al. [START_REF] Ninin | A Reliable Affine Relaxation Method for Global Optimization[END_REF]. They have proposed an interval Branch and Bound, called here IBBA+, that uses constraint propagation and a sophisticated variant of affine arithmetic. From their benchmark of 74 polynomial and non polynomial systems (without trigonometric operators), we have extracted the 27 ones that required more than 1 second to be solved by the simplest version of IbexOpt (column 4). In the extended paper [START_REF] Araya | A Contractor Based on Convex Interval Taylor[END_REF], a table shows the 11 systems solved by this first version in a time comprised between 1 and 11 seconds. Table 1 includes the 13 systems solved in more than 11 seconds. 5 Three systems (ex6 2 5, ex6 2 7 and ex6 2 13) are removed from the benchmark because they are not solved by any solver. The reported results have been obtained on a same computer (Intel X86, 3Ghz).

We have implemented the different algorithms in the Interval-Based EXplorer Ibex [START_REF] Chabert | Contractor Programming[END_REF]. Reference [START_REF] Trombettoni | Inner Regions and Interval Linearizations for Global Optimization[END_REF] details how our interval B&B, called IbexOpt, handles constrained optimization problems by using recent and new algorithms. Contraction steps are achieved by the Mohc interval constraint propagation algorithm [START_REF] Araya | Exploiting Monotonicity in Interval Constraint Propagation[END_REF] (that also lower bounds the range of the objective function). The upper bounding phase uses original algorithms for extracting inner regions inside the feasible search space, i.e., zones in which all points satisfy the inequality and relaxed equality constraints. 6 The cost of any point inside an inner region may improve the upper bound. Also, at each node of the B&B, the X-Taylor algorithm is used to produce hyperplanes for each inequality constraints and the objective function. On the obtained convex polyhedron, two types of tasks can be achieved: either the lower bounding of the cost with one call to a Simplex algorithm (results reported in columns 4 to 13), or the lower bounding and the contraction of the box, with X-NewIter (i.e., 2n + 1 calls to a Simplex algorithm; results reported in column 10) or X-Newton (columns 11, 13). The bisection heuristic is a variant of Kearfott's Smear function described in [START_REF] Trombettoni | Inner Regions and Interval Linearizations for Global Optimization[END_REF].

The first two columns contain the name of the handled system and its number of variables. Each entry contains generally the CPU time in second (first line of a multi-line) and the number of branching nodes (second line). The same precision on the cost (1.e-8) and the same timeout (TO = 1 hour) have been used by IbexOpt and IBBA+. 7 Cases of memory overflow (MO) sometimes occur. For each method m, the last line includes an average gain on the different systems. For a given system, the gain w.r.t. the basic method (column 4) is CP U time (Rand) CP U time(m) . The last 10 columns of Table 1 compare different variants of X-Taylor and X-Newton. The differences between variants are clearer on the most difficult instances. All use Hansen's variant to compute the interval gradient (see Section 3.2). The gain is generally small but Hansen's variant is more robust: for instance ex 7 2 3 cannot be solved with the basic interval gradient calculation.

In the column 3, the convexification operator is removed from our interval B&B, which underlines its significant benefits in practice. The column 4 corresponds to an X-Taylor performed with one corner randomly picked for every constraint. The next column (R+R) corresponds to a tighter polytope computed with two randomly chosen corners per inequality constraint. The gain is small w.r.t. Rand. The column 6 (R+op) highlights the best X-Taylor variant where a random corner is chosen along with its opposite corner. Working with more than 2 corners appeared to be counter-productive, as shown by the column 7 (RRRR) that corresponds to 4 corners randomly picked.

We have performed a very informative experiment whose results are shown in columns 8 (Best) and 9 (B+op): an exponential algorithm selects the best corner, maximizing the expression (4), among the 2 n ones. 8 The reported number of branching nodes shows that the best corner (resp. B+op) sometimes brings no additional contraction and often brings a very small one w.r.t. a random corner (resp. R+op). Therefore, the combination R+op has been kept in all the remaining variants (columns 10 to 14).

The column 10 (XIter) reports the results obtained by X-NewIter. It shows the best performance on average while being robust. In particular, it avoids the memory overflow on ex7 2 3. X-Newton, using ratio fp=20%, is generally slightly worse, although a good result is obtained on ex6 2 12 (see column 11).

The last three columns report a first comparison between AA (affine arithmetic; Ninin et al.'s AF2 variant) and our convexification methods. Since we did not encode AA in our solver due to the significant development time required, we have transformed IbexOpt into two variants Ibex' and Ibex'' very close to IBBA+: Ibex' and Ibex'' use a non incremental version of HC4 [START_REF] Benhamou | Revising Hull and Box Consistency[END_REF] that loops only once on the constraints, and a largest-first branching strategy. The upper bounding is also the same as IBBA+ one. Therefore we guess that only the convexification method differs from IBBA+: Ibex' improves the lower bound using a polytope based on a random corner and its opposite corner; Ibex'' builds the same polytope but uses X-Newton to better contract on all the dimensions. 9First, Ibex' reaches the timeout once more than IBBA+; and IBBA+ reaches the timeout once more than Ibex''. Second, the comparison in the number of branching points (the line Sum accounts only the systems that the three strategies solve within the timeout) underlines that AA contracts generally more than Ibex', but the difference is smaller with the more contracting Ibex'' (that can also solve ex7 2 3). This suggests that the job on all the variables compensates the relative lack of contraction of X-Taylor. Finally, the performances of Ibex' and Ibex'' are better than IBBA+ one, but it is probably due to the different implementations.

Experiments in constraint satisfaction

We have also tested the X-Newton contractor in constraint satisfaction, i.e., for solving well constrained systems having a finite number of solutions. These systems are generally square systems (n equations and n variables). The constraints correspond to non linear differentiable functions (some systems are polynomial, others are not). We have selected from the COPRIN benchmark 10 all the systems that can be solved by one of the tested algorithms in a time between 10 s and 1000 s: we discarded easy problems solved in less than 10 seconds, and too difficult problems that no method can solve in less than 1000 seconds. The timeout was fixed to one hour. The required precision on the solution is 10 -8 . Some of these problems are scalable. In this case, we selected the problem with the greatest size (number of variables) that can be solved by one of the tested algorithms in less than 1000 seconds.

We compared our method with the state of art algorithm for solving such problems in their original form (we did not use rewriting of constraints and did not exploit common subexpressions). We used as reference contractor our best contractor ACID(Mohc), an adaptive version of CID [START_REF] Trombettoni | Constructive Interval Disjunction[END_REF] with Mohc [START_REF] Araya | Exploiting Monotonicity in Interval Constraint Propagation[END_REF] as basic contractor, that exploits the monotonicity of constraints. We used the same bisection heuristic as in optimization experiments. Between two choice points in the search tree, we called one of the following contractors (see Table 2).

-ACID(Mohc): see column 3 (Ref), -X-NewIter: ACID(Mohc) followed by one call to Algorithm 1 (column 4, Xiter), -X-Newton: the most powerful contractor with ratio fp=20%, and ACID(Mohc)

as internal CP contractor (see Algorithm 2).

For X-Newton, we have tested 5 ways for selecting the corners (see columns 5-9):

-Rand: one random corner, -R+R: two random corners, -R+op: one random corner and its opposite, -RRRR: four random corners, -2R+op: four corners, i.e., two random corners and their two respective opposite ones.

We can observe that, as for the optimization problems, the corner selection R+op yields the lowest sum of solving times and often good results. The last line of Table 2 highlights that all the 24 systems can be solved in 1000 s by X-Newton R+op, while only 18 systems are solved in 1000 s by the reference algorithm with no convexification method. Each entry in Table 2 contains the CPU time in second (first line of a multi-line) and the number of branching nodes (second line). We have reported in the last column (Gain) the gains obtained by the best corner selection strategy R+op as the ratio w.r.t. the reference method (column 3 Ref), i.e., CP U time(R+op) CP U time (Ref ) . Note that we used the inverse gain definition compared to the one used in optimization (see 5.1) in order to manage the problems reaching the timeout. We can also observe that our new algorithm X-Newton R+op is efficient and robust: we can obtain significant gains (small values in bold) and lose never more than 39% in CPU time.

We have finally tried, for the scalable systems, to solve problems of bigger size. We could solve Katsura-30 in 4145 s, and Yamamura1-16 in 2423 s (instead of 33521 s with the reference algorithm). We can remark that, for these problems, the gain grows with the size.

Conclusion

Endowing a solver with a reliable convexification algorithm is useful in constraint satisfaction and crucial in constrained global optimization. This paper has presented the probably simplest way to produce a reliable convexification of the solution space and the objective function. X-Taylor can be encoded in 100 lines of codes and calls a standard Simplex algorithm. It rapidly computes a polyhedral convex relaxation following Hansen's recursive principle to produce the gradient and using two corners as expansion point of Taylor: a corner randomly selected and the opposite corner. This convex interval Taylor form can be used to build an eXtremal interval Newton. The X-NewIter variant contracting all the variable intervals once provides on average the best performance on constrained global optimization systems. For constraint satisfaction, both algorithms yield comparable results.

Compared to affine arithmetic, preliminary experiments suggest that our convex interval Taylor produces a looser relaxation in less CPU time. However, the additional job achieved by X-Newton can compensate this lack of filtering at a low cost, so that one can solve one additional tested system in the end. Therefore, we think that this reliable convexification method has the potential to complement affine arithmetic and Quad.

Algorithm 1 X

 1 -NewIter (f , x, [x]): [x] for j from 1 to m do polytope ← polytope ∪ {X-Taylor(fj,x,[x])} end for for i from 1 to n do /* Two calls to a Simplex algorithm: */ xi ← min xi subject to polytope xi ← max xi subject to polytope end for return [x]

Table 1 .

 1 Experimental results on difficult constrained global optimization systems

	1	2	3	4	5	6	7	8	9	10	11	12	13	14
	System	n	No Rand R+R R+op RRRR Best B+op XIter XNewt Ibex' Ibex" IBBA+
	ex2 1 7 20	TO 42.96 43.17 40.73 49.48	TO	TO	7.74 10.58	TO	TO 16.75
				20439 16492 15477 13200			1344	514			1574
	ex2 1 9 10	MO 40.09 29.27 22.29 24.54			9.07	9.53 46.58	103 154.02
				49146 30323 23232 19347 57560 26841	5760 1910 119831 100987 60007
	ex6 1 1	8	MO 20.44 19.08 17.23 22.66			31.24 38.59	TO	633	TO
				21804 17104 14933 14977 24204 15078 14852 13751		427468
	ex6 1 3 12	TO 1100	711	529	794	TO	TO 262.5	219	TO	TO	TO
				522036 2.7e+5 205940 211362			55280 33368
	ex6 2 6	3	TO	162	175	169	207			172	136	1033	583	1575
				172413 1.7e+5 163076 163967 1.7e+5 1.6e+5 140130 61969 1.7e+6 770332 9.2e+5
	ex6 2 8	3 97.10	121	119	110 134.7			78.1 59.3	284	274	458
			1.2e+5 117036 1.1e+5 97626 98897 1.2e+5 97580 61047 25168 523848 403668 2.7e+5
	ex6 2 9	4 25.20	33.0 36.7 35.82 44.68			42.34 43.74	455	513	523
			27892 27892 27826 27453 27457 27881 27457 27152 21490 840878 684302 2.0e+5
	ex6 2 10 6	TO 3221 2849 1924 2905			2218 2697	TO	TO	TO
				1.6e+6 1.2e+6 820902 894893 1.1e+6 8.2e+5 818833 656360
	ex6 2 11 3 10.57 19.31 7.51	7.96 10.82			13.26 11.08 41.21 11.80 140.51
			17852 24397 8498 8851 10049 5606 27016 12253 6797 93427 21754 83487
	ex6 2 12 4 2120	232	160 118.6	155			51.31 22.20	122	187 112.58
			2e+6 198156 1.1e+5 86725 90414 1.9e+5 86729 31646 7954 321468 316675 58231
	ex7 3 5 13	TO	44.7 54.9	60.3 75.63			29.88 28.91	TO	TO	TO
				45784 44443 50544 43181 45352 42453	6071 5519
	ex14 1 7 10	TO	433	445	406	489			786	938	TO	TO	TO
				223673 1.7e+5 156834 125121 1.7+5 1.1+5 179060 139111
	ex14 2 7 6 93.10 94.16 102.2	83.6 113.7			66.39 97.36	TO	TO	TO
			35517 25802 21060 16657 15412 20273 18126 12555 9723
	Sum			5564 4752 3525 5026			3767 4311	1982 1672	2963
				3.1e+6 2.2e+6 1.7e+6 1.7e+6			1.4e+6 983634 3.6e+6 2.3e+6 1.6e+6
	Gain			1 1.21	1.39	1.07			2.23	1.78
	ex7 2 3	8	MO	MO	MO	MO	MO			544	691	TO	719	TO
										611438 588791		681992

Table 2 .

 2 Experimental results on difficult constraint satisfaction problems: the best results and the gains (< 1) appear in bold

	1	2	3	4	5	6	7	8	9	10
	System	n

The aim is not to produce a convex polyhedral relaxation (which is not mentioned), but to use as expansion point the farthest point in the domain from a current point followed by the algorithm. The contraction is not obtained by calls to an LP solver but by the general purpose Gauss-Seidel without taking advantage of the convexity.

An equation f (x) = 0 can be viewed as two inequality constraints: 0 ≤ f (x) ≤ 0.

We cannot prove anything on more complicated, e.g., transcendental, functions that make the problem undecidable.

In an instance of BALANCED-3SAT, each Boolean variable xi occurs ni times in a negative literal and ni times in a positive literal. We know that BALANCED-3SAT is NP-complete thanks to the dichotomy theorem by Thomas J. Schaefer[START_REF] Schaefer | The Complexity of Satisfiability Problems[END_REF].

Note that most of these systems are also difficult for the non reliable state-of-theart global optimizer Baron[START_REF] Tawarmalani | A Polyhedral Branch-and-Cut Approach to Global Optimization[END_REF], i.e., they are solved in a time comprised between 1 second and more than 1000 seconds (time out).

[START_REF] Benhamou | Revising Hull and Box Consistency[END_REF] An equation hj(x) = 0 is relaxed by two inequality constraints: -ǫ ≤ hj(x) ≤

+ǫ.[START_REF] Bliek | Computer Methods for Design Automation[END_REF] The results obtained by IBBA+ on a similar computer are taken from[START_REF] Ninin | A Reliable Affine Relaxation Method for Global Optimization[END_REF].

We could not thus compute the number of branching nodes of systems with more than 12 variables because they reached the timeout.

We have removed the call to Mohc inside the X-Newton loop (i.e., CP-contractor=⊥) because this constraint propagation algorithm is not a convexification method.

http://www-sop.inria.fr/coprin/logiciels/ALIAS/Benches/benches.html

Acknowledgment

We would like to particularly thank G. Chabert for useful discussions about existing interval analysis results.