V Thanh Ngo

H T Diep
email: diep@u-cergy.fr

Stacked Triangular XY Antiferromagnets: End of a Controversial Issue on the Phase Transition

Keywords: numbers: 75, 10, -b General theory and models of magnetic ordering ; 75, 40, Mg Numerical simulation studies

We show in this paper by using the Wang-Landau flat-histogram Monte Carlo method that the phase transition in the XY stacked triangular antiferromagnet is clearly of first-order, confirming results from latest Monte Carlo simulation and from a nonperturbative renormalization group, putting an end to a long-standing controversial issue.

I. INTRODUCTION

Effects of the frustration in spin systems have been extensively investigated during the last 30 years. Frustrated spin systems are shown to have unusual properties such as large ground state (GS) degeneracy, additional GS symmetries, successive phase transitions with complicated nature. Frustrated systems still challenge theoretical and experimental methods. For recent reviews, the reader is referred to Ref. 1. Let us confine our discussion on the nature of the phase transition in strongly frustrated spin systems. Since the nature of the phase transition depends on the symmetry and the dimension of the system, we have to examine first its GS properties. Of course, the nature of the order parameter defined according to the system symmetry determines the properties of the phase transition. One of the most studied systems is the stacked triangular antiferromagnet (STA): the antiferromagnetic (AF) interaction between nearest-neighbor (NN) spins on the triangular lattice causes a very strong frustration. It is impossible 1 to fully satisfy simultaneously the three AF bond interactions on each equilateral triangle. The GS configuration of both Heisenberg and XY models is the well-known 120-degree structure.

The phase transition in strongly frustrated spin systems is a subject of intensive investigations in the last 20 years. Theoretically, these systems are excellent testing grounds for theories and approximations. Many well-established methods such as renormalization group (RG), high-and low-temperature series expansions etc often failed to deal with these systems. Experimentally, data on different frustrated systems show a variety of possibilities: first-order or second-order transitions with unknown critical exponents etc. (see reviews in Ref. 1). The case of XY and Heisenberg spins on the STA has been intensively studied mostly since 1987. There are good recent reviews on the subject (see for example reviews by Delamotte et al 2 ). Let us briefly recall here some main historical developments and actual situation. In the XY and Heisenberg cases, different materials give different experimental results. The anomalous dimension is found negative in many materials and in most numerical simulations, the scaling relations are violated and no universality in the exponents was found in early simulations. This situation is briefly described in the following. Kawamura 3,4 has conjectured by the two-loop RG analysis in d = 3 the existence of a new universality class for frustrated magnets. Since then there have been many other perturbative calculations with contradictory results. 5,6 From 2000, there has been a number of papers by Tessier and coworkers 7,8,9 using a nonperturbative RG study of frustrated magnets for any dimension between two and four. They recovered all known perturbative results at one loop in two and four dimensions as well as for N → ∞. They determined N c (d) for all d and found N c (d = 3) = 5.1 below which the transition is of first order. However, they found the existence of a whole region in the flow diagram in which the flow is slow. As a consequence, for N = 2, 3, they found pseudo-critical exponents in good agreement with some experimental data. This allowed them to account for the nonuniversal scaling observed in XY and Heisenberg frustrated magnets. The only problem in this nonperturbative technique is that the Hamiltonian is truncated at the beginning. Giving this fact, we have to be careful about its conclusion. As will be seen in this paper, the nonperturbative results are so far well confirmed. Early MC results on XY STA have been reviewed by Loison. 10 Until 2003, all numerical simulations found ambiguous results for this model and never a clear first-order transition. A numerical breakthrough has been realized with the results of Itakura 11 who used an improved MC renormalization-group scheme to numerically investigate the renormalization group flow of the Heisenberg and XY STA and its effective Hamiltonian which is used in the field-theoretical studies. He found that the XY STA exhibits clear first-order behavior and there are no chiral fixed points of renormalization-group flow for N=2 and 3 cases. In 2004, Peles et al 12 We study again here the XY STA with high-resolution MC technique which is very efficient specially for weak first-order transition. 14 Our aim is to try to put an end to the controversy which has been lasting for 20 years. We will recall some important numerical results in the next section.

The paper is organized as follows. Section II is devoted to the description of the model and technical details of the Wang-Landau (WL) methods as applied in the present paper.

Section III shows our results. Concluding remarks are given in section IV.

II. MONTE CARLO SIMULATION: WANG-LANDAU ALGORITHM

We consider the stacking of triangular lattices in the z direction. The spins are the classical XY model of magnitude S = 1. The Hamiltonian is given by

H = J i,j S i • S j + J ′ i,k S i • S k (1)
where S i is the XY spin at the lattice site i, i,j indicates the sum over the NN spin pairs S i and S j in a xy triangular plane, while i,k indicates that of NN spin pairs between adjacent planes. J and J ′ are in-plane and inter-plane interactions, respectively. We shall suppose that J = 1 (antiferromagnetic) and J ′ = -1 (ferromagnetic) in the following.

Recently, Wang and Landau 14 proposed a Monte Carlo algorithm for classical statistical models. The algorithm uses a random walk in energy space in order to obtained an accurate estimate for the density of states g(E). This method is based on the fact that a flat energy histogram H(E) is produced if the probability for the transition to a state of energy E is proportional to g(E) -1 .

At the beginning of the simulation, the density of states (DOS) is set equal to one for all energies, g(E) = 1. In general, if E and E ′ are the energies before and after a spin is flipped, the transition probability from

E to E ′ is p(E → E ′ ) = min [g(E)/g(E ′ ), 1] . (2)
Each time an energy level E is visited, the DOS is modified by a modification factor f > 0 whether the spin flipped or not, i.e. g(E) → g(E)f . In the beginning of the random walk the modification factor f can be as large as e 1 ≃ 2.7182818. A histogram H(E) records how often a state of energy E is visited. Each time the energy histogram satisfies a certain "flatness" criterion, f is reduced according to f → √ f and H(E) is reset to zero for all energies. The reduction process of the modification factor f is repeated several times until a final value f final which close enough to one. The histogram is considered as flat if

H(E) ≥ x%. H(E) (3)
for all energies, where the flatness parameter 0% < x% < 100% controls the accuracy of the estimated g(E), with increasing accuracy as x% approaches unity. H(E) is the average histogram.

Thermodynamic quantities 14,15 can be evaluated using the canonical distribution at any temperature T by P (E, T ) = g(E) exp(-E/k B T )/Z where Z is the partition function de-

fined by Z = E g(E) exp(-E/k B T ).
In this work, we consider a energy range of interest 16,17 (E min , E max ). We divide this energy range to R subintervals, the minimum energy of each subinterval is E i min for i = 1, 2, ..., R, and maximum of the subinterval i is E i max = E i+1 min +2∆E, where ∆E can be chosen large enough for a smooth boundary between two subintervals. The Wang-Landau algorithm is used to calculate the relative DOS of each subinterval (E i min , E i max ) with the modification factor f final = exp(10 -9 ) and flatness criterion x% = 95%. We reject the suggested spin flip and do not update g(E) and the energy histogram H(E) of the current energy level E if the spin-flip trial would result in an energy outside the energy segment. The DOS of the whole range is obtained by joining the DOS of each subinterval (E i min + ∆E, E i max -∆E). appears. This is a sufficient condition, not a necessary condition, for a first-order transition.

III. RESULTS

We

We give here the values of T c for a few sizes: T c = 1.458270, 1.457878, 1.457642, 1.457537

for N = 48, 84, 96, 120, respectively. Note that this result is in excellent agreement with earlier MC simulations 11,12,13 using less sophisticated methods. To explain why standard MC methods without histogram monitoring (see for example Ref.

3) fail to see the first order character, let us show in Fig. 3 the energy vs T obtained by averaging over states obtained by the WL method for N = 120. We see here that even at this big size, the average energy does not show a discontinuity as in a strong first-order transition: the averaging over all states erases away the bimodal distribution seen in the energy histogram at the transition temperature. Therefore care should be taken to avoid such problems due to averaging in MC simulations. We note that the distance between to peaks in Fig. 2, i. e. the latent heat, is ≃ 0.009 in agreement with earlier works.

IV. CONCLUDING REMARKS

We have studied in this paper the phase transition in the XY STA by using the flat histogram technique invented by Wang and Landau. The method is very efficient because it helps to overcome extremely long transition time between energy valleys in systems with a possible first-order phase transition. We found that the transition is clearly of first-order confirming therefore recent MC results using less efficient techniques. These results put

  have used a continuous model to study the XY STA by MC simulation. They found evidence of a first-order transition. In 2006, Kanki et al 13 , using a microcanonical MC method, have found a first-order signature of the XY STA. While these recent simulations have demontrated evidence of first-order transition for the XY STA in agreement with the nonperturbative RG analysis, all of them suffer one or two uncertain aspects: the work of Itakura has used a truncated Hamiltonian, the work of Peles et al has used standard MC methods and the work of Kanki et al used a traditional microcanonical MC technique. At present, we have a very high-performance technique at hand for weak first-order transitions. This is a very good opportunity to test it on the XY STA and to say a last word on the nature of the phase transition of this system by using the full Hamiltonian, confirming or rejecting the nonperturbative RG and recent MC results. That is the purpose of this work.

  Figs.1 and 2, respectively. As seen, for N = 48, the peak, though very large, does not show yet a double-maximum structure. Only from N = 90 that the double-peak structure clearly

FIG. 2 :

 2 FIG. 2: Energy histograms for N = 84, 96, 120 at T c indicated on the figure.

  Energy histograms for N = 48 at T c indicated on the figure.

	( ) P E
	1				T = 1.45827
	0.8
	0.6
	0.4
	0.2
	0 -1.2	-1.19	-1.18	-1.17	-1.16	E	-1.15
	FIG. 1:
				11,12,13


K Akabli

H T Diep
email: diep@u-cergy.fr

Effects of Ferromagnetic Magnetic Ordering and Phase Transition on the Resistivity of Spin Current

Keywords:

It has been shown experimentally a long time ago that the magnetic ordering causes an anomalous behavior of the electron resistivity in ferromagnetic crystals.

Phenomenological explanations based on the interaction between itinerant electron spins and lattice spins have been suggested to explain these observations. We show by extensive Monte Carlo simulation that this behavior is also observed for the resistivity of the spin current calculated as a function of temperature (T ) from low-T ordered phase to high-T paramagnetic phase in a ferromagnet. We show in particular that across the critical region, the spin resistivity undergoes a huge peak. The origin of this peak is shown to stem from the formation of magnetic domains near the phase transition. The behavior of the resistivity obtained here is compared to experiments and theories. A good agreement is observed.

I. INTRODUCTION

The interplay between magnetic properties and transport behavior has been a subject of continuous experimental investigations. 1,2,3 It has been shown in these works that the resistivity has an anomalous peak in the magnetic transition temperature region. De Gennes and Friedel, 4 and later Fisher and Langer, 5 have suggested that the peak at the transition is due to the magnetic short-range interaction. Various analytical methods such as mean-field approximations 6,7 have been employed to explain experimental observations. The problem of spin-dependent transport has been also extensively studied in magnetic thin films and multilayers. The so-called giant magnetoresistance (GMR) was discovered experimentally twenty years ago. 8,9 Since then, intensive investigations, both experimentally and theoretically, have been carried out. 10,[START_REF]Solid State Physics[END_REF] Experimental observations show that when the spin of an itinerant electron is parallel to the spins of the environment it will go through easily while it will be stopped if it encounters an antiparallel spin medium. The resistance is stronger in the latter case resulting in a GMR. Although many theoretical investigations have been carried out, to date very few Monte Carlo (MC) simulations have been performed regarding the temperature dependence of the dynamics of spins participating in the current. In a recent work, [START_REF] Akabli | [END_REF] we have investigated by MC simulations the effects of magnetic ordering on the spin current in magnetic multilayers. Our results are in agreement with measurements. 13 We study in this paper the transport of itinerant electrons traveling inside a bulk ferromagnetic crystal using extensive MC simulations. The electron spin is supposed to be the Ising model. Various interactions are taken into account, in particular interaction between itinerant spins and interaction between itinerant spins and lattice spins.

The paper is organized as follows. Section II is devoted to the description of our model and the rules that govern its dynamics. In section III, we describe our MC method and discuss the results we obtained. These results are in agreement with experiments and theories.

Concluding remarks are given in Section IV.

II. MODEL

We consider in this paper a bulk ferromagnetic crystal. We use the Ising model and the face-centered cubic (FCC) lattice with size 4N x × N y × N z . Periodic boundary conditions (PBC) are used in the three directions. Spins localized at FCC lattice sites are called "lattice spins" hereafter. They interact with each other through the following Hamiltonian:

H l = -J i,j S i • S j , (1)
where S i is the Ising spin at lattice site i, i,j indicates the sum over every nearest-neighbor (NN) spin pair (S i , S j ), J(> 0) being the NN interaction.

In order to study the spin transport in the above system, we consider a flow of itinerant spins interacting with each other and with the lattice spins. The interaction between itinerant spins is defined as follows,

H m = - i,j K i,j s i • s j , (2)
where s i is the Ising spin at position r i , and i,j denotes a sum over every spin pair (s i , s j ).

The interaction K i,j depends on the distance between the two spins, i.e.,

r ij = | r i -r j |.
A specific form of K i,j will be chosen below. The interaction between itinerant spins and lattice spins is given by

H r = - i,j I i,j s i • S j , (3)
where the interaction I i,j depends on the distance between the itinerant spin s i and the lattice spin i . For the sake of simplicity, we assume the same form for K i,j and I i,j , namely, K i,j = K 0 exp(-r ij ) and I i,j = I 0 exp(-r ij ), where K 0 and I 0 are constants.

The procedure used in our simulation is described as follows. First we study the thermodynamic properties of the bulk system alone, i.e., without itinerant spins, using Eq. ( 1). We perform MC simulations to determine quantities as the internal energy, the specific heat, layer magnetizations, the susceptibility, ... as functions of temperature T . [START_REF] Binder | Monte Carlo Simulation in Statistical Physics[END_REF] From these physical quantities we determine the critical temperature T c below which the system is in the ordered phase. We show in Fig. 1 the lattice magnetization versus T .

Once the lattice has been equilibrated at T , we inject N 0 itinerant spins into the system.

The itinerant spins move into the system at one end, travel in the x direction, escape the system at the other end to reenter again at the first end under PBC. Note that PBC are used to ensure that the average density of itinerant spins remains constant during the time (stationary regime). The dynamics of itinerant spins is governed by the following interactions: i) an electric field E is applied in the x direction. Its energy is given by

H E = -E•v
, where v is the velocity of the itinerant spin; ii) a chemical potential term which depends on the concentration of itinerant spins within a sphere of radius D 2 ("concentration gradient" effect). Its form is given by H c = Dn(r), where n(r) is the concentration of itinerant spins in a sphere of radius D 2 centered at r. D is a constant taken equal to K 0 for simplicity; iii) interactions between a given itinerant spin and lattice spins inside a sphere of radius D 1 (Eq. 3); iv) interactions between a given itinerant spin and other itinerant spins inside a sphere of radius D 2 (Eq. 2).

Let us consider the case without an applied magnetic field. The simulation is carried out as follows: at a given T we calculate the energy of an itinerant spin by taking into account all the interactions described above. Then we tentatively move the spin under consideration to a new position with a step of length v 0 in an arbitrary direction. Note that this move is immediately rejected if the new position is inside a sphere of radius r 0 centered at a lattice spin or an itinerant spin. This excluded space emulates the Pauli exclusion principle in the one hand, and the interaction with lattice phonons on the other hand. If the new position does not lie in a forbidden region of space, then the move is accepted with a probability given by the standard Metropolis algorithm. [START_REF] Binder | Monte Carlo Simulation in Statistical Physics[END_REF]

III. MONTE CARLO RESULTS

We let N 0 itinerant spins travel through the system several thousands times until a MC steps per spin and we compute statistical averages over 10 6 MC steps per spin. Taking J = 1, we obtain T c ≃ 9.75 for the estimate of the critical temperature of the lattice spins (see Fig. 1).

In Fig. 2 we sketch the characteristic traveling length λ computed after a fixed lapse of time as a function of temperature T . As can be seen, λ is very large for T < T c . We note that there is a small depression in the transition region. We will show below that this has an important consequence on the spin current: the resistance undergoes a cusp. We define the resistance R as R = 1 n , where n is the number of itinerant spins crossing a unit area perpendicular to the x direction per unit of time. In Fig. 3 we show the resistance R as a function of temperature.

There are several striking points observed here. First, R is very low in the ordered phase and large in the paramagnetic phase. Below the transition temperature, there exists a single large cluster with some isolated defects, so that any itinerant spin having the parallel orientation goes through the lattice without hindrance. The resistance is thus very small. Secondly, R exhibits a cusp at the transition temperature, the existence of which is due to the critical fluctuations in the phase transition region. De Gennes and Friedel 4 and later Fisher and Langer 5 and then Kataoka 6 among others have shown that the resistance is related to the spin correlation < S i S j > which is nothing but the magnetic susceptibility.

The latter diverges at the transition. This explains at least partially the cusp of R observed here. Another picture to explain the cusp is that when T c is approached large clusters of up (resp. down) spins form in the critical region above T c . As a result, the resistance is much larger than in the ordered phase: itinerant electrons have to steer around large clusters in order to go through the entire lattice. Thermal fluctuations are not large enough to allow the itinerant spin to overcome the energy barrier created by the opposite orientation of the clusters in this temperature region. Of course, far above T c , most clusters have a small size, so that the resistivity is still quite large with respect to the low-T phase. However, R decreases as T is increased because thermal fluctuations are more and more stronger to help the itinerant spin overcome energy barriers. Note in addition that the cluster size is now comparable with the radius D 1 of the interaction sphere, which in turns reduces the height of potential energy barriers. We have checked this interpretation by first creating an artificial structure of alternate clusters of opposite spins and then injecting itinerant spins into the structure. We observed that itinerant spins do advance indeed more slowly than in the completely disordered phase (high-T paramagnetic phase). We have next calculated directly the cluster-size distribution as a function of T (not shown) using the Kopelman algorithm. [START_REF] Hoshen | [END_REF] The result confirms the effect of clusters on the spin conductivity. The reader is referred to our previous work [START_REF] Akabli | [END_REF] for results of a multilayer case. Note that in Fig. 3 we also show the results from our theory using the Boltzmann equation in the relaxation-time approximation. The agreement is remarkable. We will show the details of this theory elsewhere. 16

IV. CONCLUDING REMARKS

We have shown in this paper first results of MC simulations on the transport of itinerant spins interacting with localized lattice spins in a ferromagnetic FCC crystal. Various interactions have been taken into account. We found that the spin current is strongly dependent on the lattice spin ordering: at low T itinerant spins whose direction is parallel (antiparallel) to the lattice spins yield a strong (weak) current. At high temperatures, the lattice spins are disordered, the current of itinerant spins is very weak and does not depend on the input orientation of itinerant spins. As a consequence, the resistance is very high at high T . We found in the transition region between low-T and high-T phases a peak of the resistance which is due to the existence of domains of lattice spins. Experiments and theories support our finding.

FIG. 1 :

 1 FIG. 1: Lattice magnetization versus temperature T . T c is ≃ 9.75 in unit of J = 1.

  steady state is reached. The parameters we use in most calculations, except otherwise stated, s = S = 1 and N x = N y = N z = 20. Other parameters are D 1 = D 2 = 1 (in unit of the FCC cell length), K 0 = I 0 = 2, N 0 = 20 3 (namely one itinerant spin per FCC unit cell), v 0 = 1, r 0 = 0.05. At each T the equilibration time for the lattice spins lies around 10 6

FIG. 2 :

 2 FIG. 2: Characteristic traveling length λ in unit of the FCC cell length versus temperature T , for 100 MC steps.

The authors are grateful to Sylvain Reynal for discussions.


R Leone

L P Lévy

P Lafarge

CPP-metrology-0 The Cooper Pair Pump as a Quantized Current Source

Keywords: numbers: 85.25.Cp, 03.65.Vf, 74.50.+r, 74.78.Na

A new charge quantization in a phase-polarized Cooper Pair Pump (CPP) is proposed, based on the topological properties of its Hamiltonian ground state over a three-dimensional parameter space P. The charge is quantized using a set of path in P covering the surface of a torus, and is a multiple of the integer Chern index c1 of this surface. This quantization is asymptotic but the pumped charge converges rapidly to the quantized value with the increase in the path frequency. The topological nature of the current makes this CPP implementation an excellent candidate for a metrological current standard.

Topological defects occur in a number of different settings. In superfluids or magnets, the ground state is described by an order parameter. Topological defects (vortices, skyrmions...) are singular points of the order parameter field [1]. A topological charge is assigned to the mapping from physical space to order parameter space. In quantum systems with few degrees of freedom, the energy levels and eigenstates may depend on a set of continuous classical parameters (or band indices). When the number of parameters is sufficient (three for a complex Hamiltonian), isolated degeneracies between two neighboring levels (or bands) can occur. The presence of a degeneracy has many physical consequences. For quantum wavefunctions, the degeneracies can then be viewed as a defect in their phase field [START_REF] Berry | Proc. R. Soc. Lond. A[END_REF]. Some classical examples are well known in polyatomic molecules [START_REF] Herzberg | [END_REF], where the nuclear coordinates form a set of semiclassical parameters (within the Born-Oppenheimer approximation) for the electronic structure. Energy manifolds can have conical intersections at isolated values of the nuclear coordinates [4,5]. When a topological charge is assigned to these degeneracies, the change of multiplicity of rotation-vibrations levels as a function of nuclear coordinates are easily understood in topological terms. In molecular magnets, isolated degeneracies have also been found for specific directions and values of the magnetic field [6,7]. The integer quantum Hall effect has also been interpreted in term of Chern indices [8,9] which are the total topological charge included within a filled Landau band. In all these examples, a topological charge is associated to an isolated degeneracy between neighboring bands, where the defect in the wavefunctions phase field occur. Its physics can be described by a "monopole" placed at the degenerate point in parameter space. The flux of the magnetic field produced by the monopole through any surface enclosing the defect is proportional to the Chern index for this surface. For a the two dimensional electron gas, the Hall conductivity of a band is σ xy = e 2 h × c 1 where c 1 is the Chern index assigned to the band [8]. The phase change of a wavefunction on a closed path is proportional to the integral of the potential vector, which is also the flux through the surface bounded by the path. This phase change, better known as Berry's phase, takes a particular value (π) for some of the physically relevant paths enclosing degeneracies. In molecular magnets, this phase change leads to destructive interferences which quench quantum tunnelling at the degeneracy points.

In this letter, the relevance of this physics to superconducting circuits is illustrated around an example, the Cooper Pair Pump (CPP), which lowest two energy levels form two bands which depend on continuous parameters (gate charges or voltages and magnetic flux). The parameter space is in this case three dimensional and degeneracies (defects) occur at isolated points in parameter space. We demonstrate how a charge quantization, proportional to the defect topological charge, can be exploited to produce a current source with metrological accuracy.

The CPP-circuit represented in Fig. 1 includes a small inductance in serie with the Cooper pair pump. The CPP has three Josephson junctions in series which define two superconducting islands. As long as the device is dominated by the Coulomb charging energy

E C = (2e) 2
4C of the islands, the topological properties of the circuit do not

E J , C E J , C E J , C ϕ = 2π(Φ/Φ 0 ) V g1 C g1 V g2 C g2 L n 1 n 2 Φ I pump
FIG. 1: The Cooper Pair Pump. depend on the precise value of the Josephson couplings E J or capacitances C, which can be taken as identical.

Let n 1 and n 2 be the number of Cooper pairs in excess on each island and

n g1 = Cg1Vg1 2e , n g2 = Cg2Vg2 2e
, the polarization charge induced by two gate voltages V g1 and V g2 . The electrostatic charges on each island are n 1n g1 and n 2n g2 , where n g1 and n g2 are classical variables tuning the islands electrostatic energies. The inductance L serves two purposes: the phase φ across the CPP is controlled by the magnetic flux Φ it threads, and it is also used as the input for the readout circuit measuring the circulating current through the pump. For a small inductance, the magnetic energy 1 2L φ -2π Φ Φ0 2 defines a sufficiently deep potential minimum: quantum fluctuations of φ across the CPP are quenched and the phase bias is set by the "classical" phase ϕ = 2π Φ Φ0 . The circuit energy can be tuned through three independent parameters n g1 , n g2 and ϕ defining the parameter space P.

Charging H C and Josephson H J Hamiltonian contribute to the CPP energy. It is convenient to write the charging energy in term of ns = n1 + n2 , the total charge on the CPP and nd = n1 -n2 , the charge imbalance between the two islands. Their canonical conjugates, the phase variables Θs and Θd enter in the Josephson Hamiltonian. In term of ns and nd the charging Hamiltonian is

H C = E C ns -n gs 2 + κ 0 nd -n gd 2 , (1)
where

E C = (2e) 2 4C
, n gs = n g1 + n g2 , n gd = n g1n g2 are the corresponding polarization charges and κ 0 = 1 3 is a capacitance ratio. Each eigenstate |n s , n d is the ground state of H C inside an hexagonal area h (ns,n d ) , centered at the point (n gs = n s , n gd = n d ), in the n gs -n gd plane (shown as the base plane in Fig. 3). Its boundaries are electrostatic degeneracy lines between "neighboring" charge states, and the vertices are points of triple degeneracies. The vertices {T 1 } = (n gs = 2 3 , n gd = 0) and {T 2 } = (n gs = 4 3 , n gd = 0) form the unit cell of this hexagonal lattice. With an appropriate gauge choice, the Josephson Hamiltonian H J depends on the phase bias ϕ as

H J = -2E J cos Θs cos Θd -E J cos(2 Θd + ϕ) .
(

) 2
The total Hamiltonian is 2π-periodic in ϕ. Lattice translations in the hexagonal lattice (n gs -n gd ) connect equivalent but physically different states. Because the parameter space P is three dimensional, H J lifts the electrostatic degeneracies except at isolated points where the ground level is still degenerate, in P (a consequence of the von Neumann-Wigner theorem [10]). These points occur at the vertices {T 1 } and {T 2 } of the hexagonal lattice for ϕ = π, as illustrated in Fig. 2. In realistic implementation of the Cooper pair pump, it is favorable to keep the ratio between Josephson and charging energy β = EJ EC small, an assumption made throughout this paper.

In the limit of zero Josephson coupling, the electrostatic states |0, 0 , |1, 1 and |1, -1 are degenerate at the vertex {T 1 }. For finite Josephson couplings, the ground state remains degenerate at this point when ϕ = π. The two degenerate ground states at T 1 = {T 1 , π} are

|g ± = 1 √ 3 |0, 0 ∓ √ 3 ∓ 1 2 |1, 1 ± √ 3 ± 1 2 |1, -1 (3)
with energy E g (T 1 ) = 4 9 E C -1 2 E J while the first excited state |e is well separated from the ground state (E e = E g + 3EJ

2 ). For a small deviation δR = σ = n gs - 2 3 , δ = n gd , ψ = ϕπ from degeneracy T 1 , a two-level approximation (|g ± ) is appropriate for small β. The eigenenergies

E ± (δR) = E g ± 2 3 E C σ 2 + δ 2 3 + 3 16 β 2 ψ 2 display a characteristic conical intersection at T 1 . After rescaling the deviation δR from T 1 , as b x = 4 3 E C σ, b y = -1 √ 3 E J ψ, b z = 4 4 √
3 δ, the twolevel Hamiltonian maps onto an isotropic spin- 1 2 problem,

Ĥ±(T1) (R) = E g 1 + 1 2 σ • b(R). (4)
This makes the conical nature of the intersection explicit (E ± (δR) = E g ± |b| 2 ), and spinor eigenstates can be specified using the direction of b rather than δR.

The topological nature of this degeneracy becomes clear when considering the adiabatic evolution of quantum states along different closed paths Γ "around" this point. One of them is represented in Fig. 3 as the time evolution of a point R(t) in P. If |g ± (R(t)) are instantaneous eigenstates of the Hamiltonian, after a period T around any closed path |g ± (T ) = e -i(η±-γ±) |g ± (0) depends on two phases, η ± = 1 T 0 E g± (R(t))dt, the dynamical phase and γ ± = i Γ A ± (R)•dR, the geometrical phase. It is expressed as an integral of the vector potential A ± (R) = g ± (R)|∇g ± (R) over the path. If the closed paths considered are sufficiently close to the degeneracy T 1 where the two-level approximation is valid, Berry's geometrical phase can also be computed using the spin- 1 2 representation as γ ± = i Γ A ± (b) • db, where the vector potential is azimuthal in this textbook example:

A ± (b) = ± cos θ∓1
2|b| sin θ e φ (the Euler angles θ, φ specify the direction of b). Its line integral over a closed path Γ is γ ± (Γ) = ∓ Ω(Γ) 2 , where Ω(Γ) is the solid angle seen from the degeneracy. The vector potentials A ± (b) can be recognized as the gauge field of Dirac monopoles of strength ∓ 1 2 placed at T 1 . They produce a radial magnetic field B ± = ∓ b 2|b| 3 which flux is ∓2π when integrated over any surface enclosing this single degeneracy. When normalized to 2π, this flux gives the Chern index c (±) 1 = ∓1 of the surface with respect to the eigenstates |g ± and one assignes to the degeneracy T 1 topological charges q ± (T 1 ) = ∓1. This topological picture is not restricted to the immediate neighborhood of T 1 where a two-level approximation can be made. Let us consider the cylindrical surface represented in Fig. 3. Since the top and bottom planes ϕ = 0 and ϕ = 2π are physically equivalent, this cylinder has the topology of a torus, and its surface is the constant radius |R| = ρ cylinder. The magnetic flux of B -= ∇ × A -through this closed surface defines its Chern index c (-) 1 of the ground state |g -(R) [11]. This index is non-zero only if the topology is non trivial, i.e. if the surface encloses degeneracies. Here the "cylindrical surface" encloses only the isolated degeneracy T 1 and Gauss theorem guarantees that its index is identical to the one of the small sphere around T 1 just considered. For a cylinder around T 2 = {T 2 , π} the Chern indices have opposite sign and q ± (T 2 ) = ±1. For a larger surface, the Chern index is the sum of the topological charges it encloses and when no degeneracies are present, it is zero.

In a CPP, pumping of charges through the circuit is achieved when moving the ground state |g -(R) adiabatically along a one-dimensional path in P. In the literature, the paths considered so far have been closed trajectories around {T 1 } or {T 2 } at constant ϕ for which the charge transferred is not quantized. When the path is chosen to cover a surface enclosing one of the singularities, the charge transferred becomes an integer multiple of the ground state Chern index. This quantization is not only observable but can be exploited to make a current standard of metrological accuracy, in a similar way as the quantum Hall effect is a standard of resistance based on the Chern indices of Landau bands. Let us consider a helical path Γ ϕ0 on the surface of the cylinder (Fig. 3) starting from the plane ϕ = ϕ 0 and ending in the plane ϕ = ϕ 0 + 2π while making an integer number p of windings (in Fig. 3, ϕ 0 = 0 and p = 2). The pumped charge Q pump on this path has a dynamical Q dyn and a geometrical contribution Q geo . The expressions

Q dyn = 2e dη - dϕ 0 (Γ ϕ0 ), Q geo = -2e dγ - dϕ 0 (Γ ϕ0 ), (5)
proved elsewere [12], generalize earlier results [13,14] to arbitrary three dimensional closed paths in P. We first show that the geometrical contribution is quantized when averaged over the initial phase ϕ 0 of the helix. Taking two helices Γ ϕ0 and Γ ϕ0+δϕ shifted in ϕ by an infinitesimal increment δϕ, a closed path Σ(ϕ 0 ) on the surface of the cylinder (shown in Fig. 3) can be constructed by adding two small opposite "vertical" segments connecting the paths Γ ϕ0 and Γ -1 ϕ0+δϕ which line integrals cancel. Using Eq. 5, the transferred charge is

Q geo = 2e γ -(Γ ϕ0 ) + γ -(Γ -1 ϕ0+δϕ ) δϕ = 2e γ -Σ(ϕ 0 ) δϕ . (6)
Using Stokes theorem, the above numerator is also

• δS B -(R) • n dS,
where δS is the surface of the strip between the two helices (shown in Fig. 3 for ϕ 0 = 0) and B -(R) is the "magnetic" field produced by the monopole at the degeneracy. When the initial angle ϕ 0 is integrated from 0 to 2π, adding all the elementary strips δS gener- ates p times the total surface S (Fig. 3). Hence

Q geo ϕ0 = 2e 2π S B -• n dS = 2e p c (-) 1 (S) = 2e p (7)
is quantized by the Chern index c 1 of the cylindrical surface S (torus). It is not necessary to average over ϕ 0 to reach quantization. If the number of windings p in Γ is large, the path integral of A -over Γ reaches the surface integral of B -over the cylinder and Q geo p→∞ = 2e p. The increase of accuracy with p being roughly an order of magnitude per additional winding, this is not only practical but also quite accurate. Using the same procedures, the dynamical charge Q dyn averages out to zero because η -(Γ ϕ0 ) is a 2π-periodic function of ϕ 0 since the Hamiltonian is 2π-periodic in ϕ. Its magnitude

≃ T E 3 J E 2 C
depends critically on the ratio β = EJ EC , however the periodic and odd-dependence of Q dyn with ϕ 0 guarantees that its ϕ 0 -average is zero. In realistic situations, ϕ 0 may not be completely controlled, or the path may not close perfectly because of flux-noise. It is therefore instructive to plot the evolution of the peak amplitude of Q dyn as a function of the number of winding p in a ϕ-period. In Fig. 4, this peak amplitude looses four order of magnitudes when p is increased by 10, using a ratio β = 0.5 an order of magnitude larger than its optimal value. For smaller β, the decrease is even more dramatic. Hence, the CPP delivers a very accurate DC current

I pump = Q pump ν p→∞ = 2e ν , ( 8
)
where ν is the winding frequency. We now consider the possible sources of errors in a practical implementation of this circuit as a current standard. Low frequency noise in the parameters will introduce some jitter in the path covered in P: the surface generated by the path will not be a perfect cylinder. On the other hand, as long as the jitter amplitude is small compared to the radius, the surface generated still encloses the topological charge, and no significant error in quantization should result. The most detrimental source of error comes from inelastic transitions between the ground state and the first excited state in the adiabatic evolution, because their topological charges are opposite. This means that if the time spent in the first exited state is τ in a period, the relative error in the transferred charge will be -4e τ T . There are two ways this can happen. For each winding around T 1 , three saddle points between the ground and first excited manifolds have to be crossed (along the hexagonal lines shown in Fig. 3), where the probability of Landau-Zener tunnelling is largest. For the optimal radius (≃ 1

3 ) this tunnelling probability is of order P ≃ e -( 3π2 )

2 E J E C E J
hν . To stay below the part per million level requires that the winding frequency ν to stay below 100 MHz for a Josephson coupling of 0.05E C , limiting the pumped current amplitude 2eν to 32 pA. The other source of inelastic transitions to the first excited state is the back-action of the measuring device on the CPP, which needs to be rigourously controlled.

The physics discussed in this paper is by no means limited to the CPP: as long as a circuit has three tunable parameters or more, degeneracies can occur. In this paper, we emphasized the topological quantization and its strong robustness to adiabatic parameters fluctuations, a key point for metrological applications. The realization of this metrological source is under way in our group: it is a beautiful scientific challenge for fundamental quantum electronics of real practical value. R. Leone is supported by a fellowship of the Rhône-Alpes Micro-Nano cluster.

FIG. 3 :

 3 FIG. 3: Representation of the degeneracy point {T1, ϕ = π} in the parameter space P. The cylindrical surface which encloses this point has the topology of a torus. The Chern index is defined as the flux through this surface normalized by 2π. Several paths on this cylinder are considered in the text.

FIG. 4 :

 4 FIG.4: Plot of the peak amplitude of Q dyn (Γϕ 0 ) as a function of the number of windings p in the interval [ϕ0; ϕ0 + 2π] for a charging energy of 10 meV, a Josephson coupling EJ = 0.5EC and a period of 10 -6 sec.

The authors are grateful to F. Faure for enlightening discussions.


Florent Bréhard

A Newton-like Validation Method for Chebyshev Approximate Solutions of Linear Ordinary Differential Systems

We provide a new framework for a posteriori validation of vectorvalued problems with componentwise tight error enclosures, and use it to design a symbolic-numeric Newton-like validation algorithm for Chebyshev approximate solutions of coupled systems of linear ordinary differential equations. More precisely, given a coupled differential system with polynomial coefficients over a compact interval (or continuous coefficients rigorously approximated by polynomials) and componentwise polynomial approximate solutions in Chebyshev basis, the algorithm outputs componentwise rigorous upper bounds for the approximation errors, with respect to the uniform norm over the interval under consideration.

A complexity analysis shows that the number of arithmetic operations needed by this algorithm (in floating-point or interval arithmetics) is proportional to the approximation degree when the differential equation is considered fixed. Finally, we illustrate the efficiency of this fully automated validation method on an example of a coupled Airy-like system.

Introduction

Notations. Let p be a positive integer for the ambient space R p , whose canonical basis is denoted by (e 1 , . . . , e p ). For a ring A, M p (A) denotes the set of order p square matrices, with 1 and 0 the identity and zero matrices.

The order over R is componentwise extended to a (partial) order over R p and M p (R): for all u, v ∈ R p (resp. A, B in M p (R)), u v if and only if u i v i for all i ∈ 1, p (resp. A B iff A ij B ij for all i, j ∈ 1, p ). The uniform norm of a function f , defined over

[-1, 1], is f ∞ = sup x∈[-1,1] |f (x)|.
General context. Numerical computing with functions is often done via polynomial approximations. For sufficiently smooth functions, defined over a compact interval, truncated Chebyshev series are well-known for their excellent approximation properties, efficient algorithms and software like Chebfun [START_REF] Boyd | Chebyshev and Fourier spectral methods[END_REF][START_REF] Trefethen | Approximation Theory and Approximation Practice[END_REF]. Broadly speaking, the principle of working with polynomial approximations instead of functions is analogous to using floating-point instead of real numbers. However, some applications such as safety-critical engineering or computer-assisted proofs in mathematics need effective and safe (rather than asymptotic) error enclosures. Problem statement and contributions. We present a symbolic-numeric a posteriori validation algorithm that provides componentwise and tight error enclosures for Chebyshev approximations to solutions of coupled linear ordinary differential equations (LODEs):

Y (r) + A r-1 (t) • Y (r-1) + • • • + A 1 (t) • Y + A 0 (t) • Y = G(t), (1)
of unknown Y : [-1, 1] → R p . Coefficients A i and G must be continuous functions, given as polynomials with rigorous error bounds. However, for the sake of simplicity, we mainly focus on the polynomial case, and refer to the solutions as vector-valued D-finite functions. Although such functions can be seen as vectors of (scalar) D-finite functions, the decoupling of the system followed by a possible desingularization step may produce hard to validate scalar LODEs (see Section 4). Moreover, in the nonpolynomial case, such techniques do not apply. Using an appropriate integral transform of the linear differential system, we obtain a Volterra integral equation of the second kind with polynomial kernel, whence the following problem statement: Problem 1. For a given integral equation of unknown Φ : [-1, 1] → R p :

Φ(t) + t -1 K(t, s) • Φ(s)ds = Ψ(t),
with a p-dimensional polynomial kernel K(t, s) ∈ M p (R[t, s]) and Ψ ∈ R[t] p , assuming we are given for each component Φ i of the exact solution Φ a polynomial approximation Φ • i in Chebyshev basis, compute componentwise error bounds ε i , as tight as desired:

Φ • i -Φ i Ч 1 ε i , for all i ∈ 1, p .
Here, • Ч 1 is a norm for absolutely summable Chebyshev series that upperbounds the • ∞ norm over [-1, 1] (see Section 3.1). Fixed-point methods are extensively used in the field of functional analysis and differential equations. They provide iterative approximation schemes, like Picard-Chebyshev which integrates nonlinear dynamical systems arising, for instance, in space flight mechanics problems [START_REF] Clenshaw | The solution of nonlinear ordinary differential equations in Chebyshev series[END_REF][START_REF] Bai | Modified Chebyshev-Picard iteration methods for solution of initial value and boundary value problems[END_REF]. They also underlie numerous validation methods for function space problems [START_REF] Kaucher | Self-validating numerics for function space problems: Computation with guarantees for differential and integral equations[END_REF][START_REF] Yamamoto | A numerical verification method for solutions of boundary value problems with local uniqueness by Banach's fixed-point theorem[END_REF].

Many fixed-point validation methods use the Banach fixed-point theorem. Given an equation x = T • x with T contracting of ratio λ ∈ (0, 1) over a complete metric space, and an approximation x to the exact solution x , it provides an error enclosure:

x -T • x 1 + λ x -x x -T • x 1 -λ . (2)
However, in the case we consider, x belongs to a product space, and the classical method consisting in endowing it with a global norm fails to produce componentwise tight error enclosures. This is particularly annoying when the components of the system are of different nature (e.g., position and speed) or magnitude.

Based on a new refinement with lower bounds for the Perov fixed-point theorem (a vector-valued generalization of Banach fixed-point principle), we propose a validation algorithm to solve Problem 1. It is a generalization of the validation method presented in [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF] to vectorial LODEs, within a new general framework for vector-valued fixed-point validation.

Theorem 1. Algorithms 1 and 3 solve together Problem 1 by providing componentwise error enclosures, as tight as desired.

(i) Algorithm 1 only depends on the integral equation (not on the provided approximation). It produces and rigorously bounds a Newton-like validation operator and requires O(p 3 N val 2 d) arithmetic operations. (ii) Algorithm 3 computes the error enclosures for the approximation and runs in linear time with respect to the maximum degree of the approximations Φ • i and the right-hand sides Ψ i . More precisely, its complexity is

O(p 2 d 2 N app + pN rhs + p 2 N val min(max(N app + d, N rhs ), N val ))
, where:

• N app = max i deg Φ • i and N rhs = max i deg Ψ i ; • d = 1 + max ij deg k ij (t, s);
• N val is a truncation index used to rigorously approximate the problem in finite dimension.

We assume a uniform complexity model, i.e., a unit cost for each arithmetic operation (+, -, ×, /, √ ), with, say, floating-point or interval operands.

The previous complexity estimates still involve a truncation index N val , which is directly related to how tight the desired error enclosures have to be. As detailed in Theorem 4, its minimal value ensuring a contracting Newtonlike operator is potentially exponential with respect to the magnitude of the coefficients of the integral equation, in the case of stiff LODEs for example. In practice however, this method works efficiently and fully automatically. An open source library implementing this validation method (and its extension to the nonpolynomial case) can be found here1 . It was also recently used for a space flight dynamics application [START_REF] Arantes Gilz | Validated Semi-Analytical Transition Matrix for Linearized Relative Spacecraft Dynam-ics via Chebyshev Polynomials[END_REF]. Previous work. In this context, applications of the Banach fixed-point theorem include early works [START_REF] Kaucher | Self-validating numerics for function space problems: Computation with guarantees for differential and integral equations[END_REF][START_REF] Yamamoto | A numerical verification method for solutions of boundary value problems with local uniqueness by Banach's fixed-point theorem[END_REF], where variations of Newton's method perform a posteriori validation in function spaces. More recent works developed techniques (e.g., radii polynomials [START_REF] Hungria | Rigorous numerics for analytic solutions of differential equations: the radii polynomial approach[END_REF]) to find a stable neighborhood of an approximation ϕ over which the Banach fixed-point theorem applies. They have the advantage of dealing with nonlinear problems (examples can be found in [START_REF] Lessard | Rigorous numerics for nonlinear differential equations using Chebyshev series[END_REF][START_REF] Van Den Berg | Rigorous numerics in dynamics[END_REF][START_REF] Hungria | Rigorous numerics for analytic solutions of differential equations: the radii polynomial approach[END_REF]). However, the above mentioned methods were not fully automated and little emphasis was put on their algorithmic aspects.

By contrast, [START_REF] Benoit | Rigorous uniform approximation of D-finite functions using Chebyshev expansions[END_REF] is a pioneer work towards effective methods for validation of approximations of D-finite functions in Chebyshev basis. At the cost of a more restricted class of functions, namely, D-finite functions, this article introduces a fully automated algorithm together with complexity estimates, based on a Picard iteration scheme. In line with this work, [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF] describes another algorithm based on a Newton-like method in an appropriate function space, which is easily extended to the case of continuous coefficients rigorously approximated by Chebyshev polynomials.

The above mentioned validation techniques are usually transposed to the vectorial case by fixing a norm over the vector-valued function space. However, this does not provide componentwise tight error enclosures. To overcome this limitation, we consider the notion of vector-valued (or generalized ) metric spaces and generalized contractions (or P-contractions) [START_REF] Kantorovich | Functional analysis in partially ordered spaces[END_REF][START_REF] Robert | Étude et utilisation de normes vectorielles en analyse numérique linéaire[END_REF][START_REF] Ortega | Iterative solution of nonlinear equations in several variables[END_REF]. The Perov fixed-point theorem [START_REF] Kantorovich | Functional analysis in partially ordered spaces[END_REF][START_REF] Perov | On the Cauchy problem for a system of ordinary differential equations[END_REF] is a natural extension of the Banach fixed-point theorem and provides componentwise upper bounds for the approximation error. Several works applied this theorem in various settings, for example [START_REF] Yamamoto | A unified derivation of several error bounds for Newton's process[END_REF] for the Newton method or [START_REF] Agarwal | Contraction and approximate contraction with an application to multi-point boundary value problems[END_REF][START_REF] Precup | The role of matrices that are convergent to zero in the study of semilinear operator systems[END_REF][START_REF] Nica-Bolojan | Fixed point methods for nonlinear differential systems with nonlocal conditions[END_REF] for ODEs with nonlocal conditions. To the best of our knowledge, however, none of these works investigate the existence of lower bounds, nor address validation problems. Outline. Section 2 introduces a general framework for componentwise fixed-point validation in generalized metric spaces. In Section 3, we design the Newton-like validation algorithm for Chebyshev approximations of vectorvalued D-finite functions. Finally, Section 4 details the validation of a twodimensional highly oscillating system. For completeness, we also provide a comparison with a decoupling technique that boils down to solving scalar LODEs.

A Framework for Vector-Valued Validation Problems

We address the general problem of componentwise validating an approximation x • to the exact solution x of a fixed-point equation x = T•x. Section 2.1 gives a rigorous definition of "several components and norms" with the notion of generalized metric spaces, leading to the Perov fixed-point theorem. Section 2.2 presents a new result that complements the Perov theorem with lower bounds on the componentwise approximation errors.

A toy example in the plane illustrates the vector-valued validation framework. Consider the trigonometric equation sin 3 ϑ + cos 3ϑ = 0 for ϑ ∈ R. By introducing c = cos x and s = sin x, this is equivalent to finding the roots of the following polynomial system in the plane (c, s):

F • (c, s) = s 3 + 4c 3 -3c c 2 + s 2 -1 = 0. (3)
Let x = (c , s ) be an exact solution and x • := (c • , s • ) = (0.84, 0.55) an approximation of it. In order to validate this solution with respect to a given norm • on R 2 , we define a Newton-like operator T•(c, s) = (c, s)-A•F•(c, s)

with A := 0.25 -0.20 -0.37 1.2 ≈ (DF x • ) -1 ∈ M 2 (R) an approximate inverse of the Fréchet derivative DF x • of F at x • .
Since A is injective, its fixed points are exactly the roots of F. In this example, F is nonlinear, so one must find a stable closed neighborhood over which T is contracting, for the Banach theorem to apply. It suffices to determine a radius r > 0 satisfying the following two conditions:

(i) λ := sup x-x • r 1 -A • DF x < 1; (ii) x • -T • x • + kr r.
If such a radius exists, then by the Banach fixed-point theorem, we have

x • -x A • F • x /(1 -λ).
However, such a bound captures a "global" error, which may not be what we expect, if, for example, the two components are of different nature (e.g., position and velocity), or differ by several orders of magnitude.

Generalized Metric Spaces and Perov Fixed-Point Theorem

Definition 1. Let X be a set (resp. E a linear space). A function d :

X × X → R p + (resp. • : E → R p +
) is a vector-valued or generalized metric (resp. norm) if for all x, y, z in X or E and λ ∈ R:

• d(x, y) = 0 iff x = y, resp. x = 0 iff x = 0; • d(x, y) = d(y, x), resp. λx = |λ| x ; • d(x, y) d(x, z) + d(z, y),
resp.

x + y x + y . Then (X, d) (resp. (E, • )) is a vector-valued or generalized metric space (resp. linear space).

A straightforward example is the product of p metric spaces (X i , d i ), i ∈ 1, p (resp. p normed linear spaces (E, • i )) and the vector-valued metric d(x, y) = (d 1 (x 1 , y 1 ), . . . , d p (x p , y p )) (resp. the vector-valued norm x = ( x 1 1 , . . . , x p p )). Remark 1. A vector-valued metric space (respectively a vector-valued normed linear space) can be trivially seen as a metric space (respectively a normed linear space) by taking the maximum of all the components of the vector-valued metric (respectively norm). We therefore recover all the useful topological notions of convergence, limit, neighborhood, completeness, etc.

In the context of vector-valued metric spaces, the notion of contracting map needs to be generalized. Let M →0 p (R) ⊆ M p (R) denote the convergent to zero matrices, that is the matrices M such that M k → 0 as k → ∞. Equivalently, these are matrices M with spectral radius ρ(M ) < 1. Then,

M →0 p (R + ) = M →0 p (R) ∩ M p (R + )
denotes those among them with nonnegative coefficients. Definition 2. Let (X, d) be a vector-valued metric space and T : X → X an operator.

• T is Λ-Lipschitz for some Λ ∈ M p (R + ) if: d(T • x, T • y) Λ • d(x, y),
for all x, y ∈ X.

• If moreover Λ is convergent to 0 (Λ ∈ M →0 p (R + )
), then T is said to be a generalized contraction.

Using these definitions, the Perov fixed-point theorem2 is a generalization of the Banach fixed-point theorem.

Theorem 2 (Perov). Let (X, d) be a complete vector-valued metric space and T : X → X a generalized contraction with a Lipschitz matrix Λ ∈ M →0 p (R + ). Then:

(i) T admits a unique fixed-point x ∈ X;

(ii) for every x • ∈ X, the iterated sequence defined by x 0 = x • and x n+1 = T • x n converges to x with the following upper bound on the approximation error:

d(x n , x ) Λ n • (1 -Λ) -1 • d(x • , T • x • ), for all n ∈ N. (4)
A proof of this theorem is given in Appendix A.1 or [START_REF] Ortega | Iterative solution of nonlinear equations in several variables[END_REF]. Perov theorem applied to the toy example. Endowing R 2 with the vector-valued norm (c, s) := (|c|, |s|) does not change the definition of T. The two conditions needed to apply the Banach fixed-point theorem are adapted to the Perov theorem as follows. Choose a multi-radius r = (r 1 , r 2 ) such that

(i) Λ := sup x-x • r |(DT x ) ij | 1 i,j 2 satisfies ρ(Λ) < 1; (ii) x • -T • x • + Λ • r r.
For r = (0.005, 0.005), one obtains:

Λ = 5.81 1.31 5.63 3.40 • 10 -2 , ρ(Λ) = 7.57 • 10 -2 ,
which satisfies (i) and (ii). Hence, Theorem 2 gives:

|c • -c | 2.90 • 10 -3 , |s • -s | 3.65 • 10 -3 .
To assess the tightness of these bounds, we provide lower bounds on the componentwise approximation errors.

Lower Bounds and Error Enclosures

Let ε = d(x • , x ) ∈ R p + be the vector of unknown errors and η = d(x • , T•x • ) ∈ R p + . By the triangle inequality, ε is circumscribed into a polytope of R p + : (1 -Λ) • ε η, (1 + Λ) • ε η, ε 0. (5)
The first inequality gives the upper bounds ε + = (1 -Λ) -1 • η, as stated by Theorem 2 (with n = 0). However, the second one does not directly give the desired lower bounds, say ε -, because the inverse

(1 + Λ) -1 = k 0 (-Λ) k is not nonnegative in general.
It is clear that each ε - i is given by the i-th coordinate of some vertex of this polytope. Instead of testing its 2 p vertices, the following theorem identifies the correct one.

Theorem 3 (Lower bounds for the Perov theorem). With the above notations, for each i ∈ 1, p , the lower bound ε - i on the i-th component ε i of the approximation error of x • to x is given by the i-th component of the vertex defined by the intersection of the i-th lower-bound constraint together with all the j-th upper-bound constraints with j = i from (5). Formally:

ε i ε - i with ε - i = e T i • (1 -D i • Λ) -1 • η,
where D i is the order p diagonal matrix defined by

(D i ) ii = -1 and (D i ) jj = 1 for j = i.
Remark 2. Contrary to the one-dimensional case, ε - i may be negative (we then round it to 0): the overestimation factor of ε + i provided by Theorem 2 is not controlled. A tighter enclosure can be obtained with a more contracting T (see Appendix A.2).

Proof of Theorem 3. Among the Inequalities (5), take the p upper-bound constraints and replace the i-th one by the corresponding lower-bound constraint. Multiply these p -1 upper-bound constraints by -1 to obtain the following system of inequalities:

(Λ -D i ) • ε -D i • η. (6)
From Lemma 1, Λ -D i is nonsingular and its inverse has nonnegative coefficients on its i-th row. Hence we can multiply ( 6) by (Λ -D i ) -1 and only keep the resulting i-th constraint:

ε i = e T i • (Λ -D i ) -1 • (Λ -D i ) • ε e T i • (Λ -D i ) -1 • (-D i ) • η = e T i • (1 -D i • Λ) -1 • η.
Lower bounds for the toy example. The polytope given by the linear constraints ( 5) is depicted in Figure 1. The top right vertex corresponds to (ε + 1 , ε + 2 ). Also, the ε - 1 (resp. ε - 2 ) is given by the top left (resp. bottom right) vertex, which is consistent with Theorem 3. This gives the following numerical enclosures:

l0pt 2.6 2.8 3 3.2 3.4 3.6 • η ε1 (×10 -3 ) ε2 (×10 -3 ) ε - 1 ε + 1 ε - 2 ε + 2 Figure 1:
Error polytope for the toy example.

ε - 1 = 2.48 • 10 -3 ε 1 = |c • -c | ε + 1 = 2.90 • 10 -3 , ε - 2 = 3.09 • 10 -3 ε 2 = |s • -s | ε + 2 = 3.65 • 10 -3 . (7)
The tightness of these enclosures is discussed in Appendix A.2. Roughly speaking, the ratio ε + i /ε - i depends not only on Λ (like in the univariate case), but also on η = d(x • , T • x • ).

Componentwise Validation of Chebyshev Approximations

We present the validation method to solve Problem 1. Section 3.1 contains reminders about Chebyshev approximation theory and LODEs. This leads to an efficient approximating procedure (Section 3.2). Section 3.3 presents Algorithms 1 and 2 to create and bound a Newton-like operator associated to a given vectorial LODE, then Algorithm 3 to compute componentwise error enclosures for any Chebyshev approximation (Φ • i ) 1 i p .

Reminders on Chebyshev Approximations

Chebyshev series and Ч 1 space. The Chebyshev family of polynomials is defined by the three-term recurrence T n+2 = 2XT n+1 -T n with initial terms T 0 = 1 and T 1 = X. They satisfy the fundamental trigonometric relation T n (cos ϑ) = cos(nϑ), from which we deduce some of their basic algebraic properties:

T n T m = 1 2 (T n+m + T |n-m| ), T n = T n+1 2(n + 1) - T n-1 2(n -1) = T n (n 2), (8)
and that

|T n (t)| 1 for x ∈ [-1, 1]. Let L 2 Ч = L 2 (1/ √ 1 -t 2 ) denote the space of real-valued measurable func- tions f over [-1, 1] such that 1 -1 f (t) 2 / √ 1 -t 2 dt < ∞. The inner product f, g = 1 -1 f (t)g(t)/ √ 1 -t 2 dt = π 0 f (cos ϑ)g(cos ϑ)dϑ,
defines a Hilbert space structure over L 2 Ч , for which the Chebyshev polynomials form a complete orthogonal system. To any continuous function f in this space we can associate its Chebyshev coefficients:

[f ] n = 1 π π 0 f (cos ϑ)dϑ, if n = 0, 2 π π 0 f (cos ϑ) cos(nϑ)dϑ, if n > 0.
Hence, the truncated Chebyshev series

f [N ] = п N • f := N n=0 [f ] n T n of f
is the orthogonal projection of f onto the finite-dimensional subspace spanned by T 0 , . . . , T N . In addition to the L 2 Ч convergence, and analogously to Fourier series, Chebyshev series have excellent approximation properties [START_REF] Boyd | Chebyshev and Fourier spectral methods[END_REF]. For example, if f is of class C r over [-1, 1] with r 1, then f [N ] uniformly converges to f in O(N -r ), and the convergence is even exponential for analytic functions.

We call Ч 1 the Banach space of continuous functions with absolutely summable Chebyshev series, with norm

f Ч 1 = n 0 |[f ] n |. Note that Ч 1
is analogous to the Wiener algebra A(T) of absolutely convergent Fourier series [12, §I.6]: for f ∈ Ч 1 , we have f Ч 1 = f (cos) A(T) . We obtain a Banach algebra structure:

f g Ч 1 f Ч 1 g Ч 1 .
Moreover, this norm is a safe overestimation of the uniform norm:

f Ч 1 sup -1 t 1 n 0 |[f ] n T n (t)| sup -1 t 1 |f (t)| = f ∞ .
Given an endomorphism F : Ч 1 → Ч 1 , the operator norm induced by the Ч 1 norm is given by:

F Ч 1 = sup f Ч 1 1 F • f Ч 1 = sup n 0 F • T n Ч 1 . ( 9
)
This corresponds to the maximum sum of the coefficients in absolute value over all columns of the matrix representation of F.

D-finite equations and integral transforms. We consider a generic pdimensional order r system of LODEs over the compact interval [-1, 1]:

Y (r) + A r-1 (t) • Y (r-1) + • • • + A 1 (t) • Y + A 0 (t) • Y = G(t), (10)
of unknown Y = (Y 1 , . . . , Y p ) : [-1, 1] → R p , with polynomial coefficients A k = (a kij ) 1 i,j p ∈ M p (R[t]) and right-hand side G = (G 1 , . . . , G p ) ∈ R[t] p .
We also fix initial conditions at -1:

Y (i) (-1) = v i , v i ∈ R p , for all i ∈ 0, r -1 . (11)
Together, ( 1) and ( 11) form an Initial Value Problem (IVP).

Several barriers arise when working directly on a differential equation ( 1): the differentiation of Chebyshev polynomials does not admit a compact formula, whence a dense linear system to solve, and, from the theoretical point of view, the space Ч 1 is not stable under differentiation. A common way to circumvent these limitations is to apply an integral transform onto the IVP problem so as to obtain an equivalent Volterra integral equation of the second kind over [-1, 1]:

Φ + K • Φ = Ψ, with K • Φ(t) = t -1 K(t, s) • Φ(s)ds, (12)
with a bivariate polynomial kernel K = (k ij ) 1 i,j p ∈ M p (R[t, s]) and righthand side Ψ = (Ψ 1 , . . . , Ψ p ) ∈ R[t] p . Depending on the integral transform, the unknown function Φ = (Φ 1 , . . . , Φ p ) : [-1, 1] → R p can be either Y or one of its derivatives. For example, [START_REF] Benoit | Rigorous uniform approximation of D-finite functions using Chebyshev expansions[END_REF] acts over Y , whereas [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF] considers the last derivative Y (r) . In any case, K : Φ → t -1 K(t, s)•Φ(s)ds is a bounded linear operator from (Ч 1 ) p to itself. We may describe it by blocks K = (K ij ) 1 i,j p , where each K ij is a one-dimensional integral operator of kernel k ij (t, s). By decomposing k ij (t, s) in Chebyshev basis with respect to s, we obtain unique polynomials b ijk (t) such that

k ij (t, s) = κ ij k=0 b ijk (t)T k (s), K ij • ϕ(t) = κ ij k=0 b ijk (t) t -1 T k (s)ϕ(s)ds.
Consequently to the multiplication and integration formulas (8), the (infinite dimensional) matrix representation of [START_REF] Olver | A fast and well-conditioned spectral method[END_REF], meaning that the nonzero entries are located on the h ij first rows (horizontal band with initial entries) and the diagonal plus the first d ij upper and lower diagonals (diagonal band with diagonal entries), with 2(a)).

K ij : Ч 1 → Ч 1 has a so-called (h ij , d ij ) almost-banded structure
h ij = max 0 k κ ij deg b ijk (t) and d ij = 1 + deg k ij (t, s) = 1 + max 0 k κ ij (k + deg b ijk (t)) (see Figure

Efficient numerical solving

The integral equation ( 12) is an infinite-dimensional linear system over the Chebyshev coefficients of the unknown function Φ. The projection method (also sometimes called Galerkin method [START_REF] Gottlieb | Numerical Analysis of Spectral Methods: Theory and Applications[END_REF]) consists in truncating for a given index N app and solving the obtained finite-dimensional linear system. In our case, this can be efficiently done by taking advantage of its sparse structure.

Define the N app -th truncation of K as

K [Napp] = (K [Napp] ij ) 1 i,j p , where п Napp • K ij • п Napp (see
K [Napp] recovers a (ph, pd) almost-banded structure, where h = max ij h ij and

d = max ij d ij (see Figure 2(d)).
Hence, solving the approximate problem:

Φ + K [Napp] • Φ = Ψ requires O(p 3 N app d 2 )
operations, using the algorithm of [START_REF] Olver | A fast and well-conditioned spectral method[END_REF] for solving almostbanded linear systems.

Validation Procedure

We extend the validation procedure of [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF] to the vectorial case. We prove the main Theorem 1 in order to solve Problem 1 in two steps: (1) a Newton-like validation operator is created and bounded by Algorithm 1. This first step is independent of the approximation degree N app . (2) The error enclosure of the given approximation is computed by Algorithm 3, following Theorems 2 and 3.

Newton-like validation operator. Following the idea of Newton's method and similar approaches, Equation ( 12) is transformed into the fixed-point equation:

T • Φ = Φ, T • Φ := Φ -A • (Φ + K • Φ -Ψ), (14)
which is equivalent to (12) as soon as A : (Ч 1 ) p → (Ч 1 ) p is injective. Moreover, T is an affine operator of linear part DT = 1 -A • (1 + K). The main challenge is to efficiently compute A and bound DT (Ч 1 ) p . This is ) -1 , for some truncation order N val . Choosing N val is a tradeoff between proving T is contracting (N val must be large enough so that K -K [N val ] (Ч 1 ) p is rigorously proved to be sufficiently small) and efficiency requirements (see [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF] for heuristics to find N val ).

(a) Almost-banded structure of K ij (b) Truncation K [Napp] ij of K ij (c) K [Napp] in B p,Napp , by blocks (p = 3) (d) K [Napp] in B p,N app is almost- banded
Once N val is fixed, Algorithm 1 first computes an approximate inverse A in floating-point (lines 1-4). Since 1 + K [N val ] is almost-banded in B p,N val , its numerical inverse can be either computed with the numerically stable algorithm of [START_REF] Olver | A fast and well-conditioned spectral method[END_REF], or approximated by a (ph , pd ) almost-banded matrix [7, Algorithm 5], in O(p 3 N val (h + d )(h + d)) operations. The operator A is defined by extending A to the whole space (Ч 1 ) p by the identity.

Second, Algorithm 1 bounds a Lipschitz matrix for T, as DT (Ч 1 ) p = ( (DT) ij Ч 1 ) 1 i,j p , block by block, using the triangle inequality:

DT Ч 1 1 -A • (1 + K [N val ] ) Ч 1 + A • (K -K [N val ] ) Ч 1 . ( 15
)
The first part of ( 15) is the approximation error, measuring how far A is from the inverse of 1 + K [N val ] . This is straightforwardly bounded as Λ A by Algorithm 1 (lines 5-9) using O(p 3 N val (h + d )(h + d)) interval arithmetic operations, and the resulting bound takes into account all sources of errors: rounding errors, sparse approximation, etc. Since only additions and multiplications of matrices are involved, the use of interval arithmetics is not critical. However, if needed, the underlying floating-point precision can be increased.

The second part of ( 15) is the truncation error, because the truncated operator K [N val ] only approximates K. Let E ij be the (i, j) block of E := A • (K -K [N val ] ):

E ij = p k=1 A ik • (K kj -K [N val ] kj ). (16)
Algorithm 1 (lines 10-16) computes Λ T E (Ч 1 ) p by blocks, with the triangle inequality: each subterm of ( 16) is rigorously bounded by Algorithm 2. This algorithm, detailed below, requires O((h

+ d )(h + d) 2 ) in- terval arithmetic operations. Hence the computation of Λ T is in O(p 3 (h + d )(h + d) 2 ).
Finally, Algorithm 1 computes Λ = Λ A + Λ T and checks that this Lipschitz matrix is convergent to zero, in which case the constructed Newton-like operator T is contracting. The eigenvalues of Λ can be safely computed with interval arithmetics, for the dimension p is usually small (typically, p 100).

Proof of Theorem 1(i).

The detailed description of Algorithm 1 above proves its correctness, and the given complexity estimates for lines 1-4, 5-9 and 10-16 sum to a global complexity of O(p 3 N val (h + d )(h + d)) operations. In the worst case, when A is dense (h + d ≈ N val ), we recover the estimate of Theorem 1(i).

Truncation error bounding. From Equation ( 16), one needs to bound

A ik •(K kj -K [N val ] kj ) Ч 1
, where A ik is the extension to Ч 1 of the order N val +1 matrix A ik by the identity if i = k, and zero otherwise. This computation is handled by Algorithm 2, which is a modification of [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF]Algorithm 6], that only treats the case i = k.

Specifically, let K denote here a one-dimensional (h, d) almost-banded integral operator, and A : Ч 1 → Ч 1 the extension of an order N val + 1 matrix A by the identity or zero. We have:

A • (K -K [N val ] ) Ч 1 = sup 0 B( ), with B( ) = A • (K -K [N val ] ) • T Ч 1 .
Indices are divided into four groups, reflecting how the initial and diagonal coefficients are impacted by the action of A:

0, N val -d , N val -d + 1, N val , N val + 1, N val + d and N val + d + 1, +∞ .
The T in the first group lie in the kernel. The second and third ones are explicitly computed, yielding bounds δ (1) and δ (2) in Algorithm 2 (lines 1-7 and 8-13). For the infinite last group, B( ) is decomposed as B I ( ) + B D ( ), the contribution of the initial and diagonal coefficients. Algorithm 2 uses the efficient bounding strategy of [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF]. First, it computes the image of T N val +d+1 for the initial and diagonal coefficients (lines 16 and 22). Then, it bounds the difference between the images of T N val +d+1 and the remaining T for > N val + d + 1 to finally deduce bounds δ (3) and δ (4) (lines 17 and 23). Error enclosures. Finally, Algorithm 3 implements the validation procedure of Theorems 2 and 3 by applying the operator T to the candidate approximation Φ • , bounding the distance of the resulting polynomial to Φ • and producing componentwise error enclosures to Φ with respect to the Ч 1 norm.

Proof of Theorem 1(ii). Algorithm 3 computes Φ

• -T • Φ • = A • (Φ • + K • Φ • -Ψ). Each P k (line 1
) is a polynomial of degree at most max(N app + d, N rhs ), and computing its Chebyshev coefficients is in O(pd 2 N app + N rhs ). Then, the computation of the coefficients of each

A ik • п N val • P k (line 3) is in O((h + d ) deg(п N val • P k )) = O((h + d ) min(max(N app + d, N rhs ), N val )).
Finally, the complexity of computing the enclosures (lines 6-7) only depends on p, and is therefore negligible. The overall complexity is:

O(p 2 d 2 N app + pN rhs + p 2 (h + d ) min(max(N app + d, N rhs ), N val )),
which gives the estimate of Theorem 1(ii) when h , d ≈ N val .

Algorithm 1 Create and bound a Newton-like operator T Require: A polynomial integral operator K = (K ij ) 1 i,j p given by the (b ijk ) 1 i,j p 0 k κ ij , and a truncation order N val . Ensure: An approximate inverse A of 1 + K [N val ] and a certified Lipschitz matrix Λ for 1 -A • (1 + K), or fail if N val not large enough.

Compute an approximate inverse matrix A.

1: M = (M ij ) 1 i,j p ← 1 + K [N val ]
, by blocks 2: M ← M in basis B p,N val 3: A ← a numerical approximate inverse of M , either dense or almostbanded.

4: A = (A ij ) 1 i,j p ← A in basis B p,N val ,

by blocks

Compute the approx error Λ A = (λ A ij ) in interval arith. 5: for i = 1 to p and j = 1 to p do 6:

C ← 1 k p A ik • M kj 7: if i = j then C ← C -1 N val +1 8: λ A ij ← C Ч 1 9: end for
Compute the trunc error Λ T = (λ T ij ) in interval arith. 10: for i = 1 to p and j = 1 to p do 11:

λ T ij ← 0 12:
for k = 1 to p do 13:

δ ← Algorithm 2 on K jk , A ik and diag := (i = k).

14:

λ T ij ← λ T ij + δ 15:
end for 16: end for Compute Λ and check if T contracting. 17:

Λ ← Λ A + Λ T 18: if ρ(Λ) < 1 then 19: return A, Λ 20: else 21:
print "Fail, Λ is not convergent to 0" 22: end if Algorithm 2 Bound the truncation error Require: A polynomial (one-dimensional) integral operator K given by the (b k ) 0 k κ , a truncation order N val , a N val + 1 order square matrix A, and a Boolean diag. Ensure: An upper bound δ for A•(K-K [N val ] ) Ч 1 , where A is the extension of A to the whole space Ч 1 by the identity if diag = true, and by zero otherwise.

All operations are performed in interval arithmetics Compute δ (1) sup ∈ N val -d+1,N val B( ) (2) then δ (2) ← P Ч 1 13: end for Compute δ (3) sup N val +d+1 B D ( )

1: δ (1) ← 0 2: if diag then 3: for = N val -d + 1 to N val do 4: P ← (1 -п N val ) • K • T 5: if P Ч 1 > δ (1) then δ (1) ← P Ч 1 6: end for 7: end if Compute δ (2) sup ∈ N val +1,N val +d B( ) 8: δ (2) ← 0 9: for = N val + 1 to N val + d do 10: P ← A • п N val • K • T 11: if diag then P ← P + (1 -п N val ) • K • T 12: if P Ч 1 > δ
14: 0 ← N val + d + 1 and B ← κ k=0 b k Ч 1 15: if diag then 16: P ← (1 -п N val ) • K • T 0 17: δ (3) ← P Ч 1 + (κ+1)B ( 0 -(κ-1)) 2 18: else 19: δ (3) ← 0 20: end if Compute δ (4) sup N val +d+1 B I ( ) 21: B ← κ k=0 A • b k Ч 1 22: P ← A • п N val • K • T 0 23: δ (4) ← P Ч 1 + (κ+1) 3 B ( 2 0 -(κ+1) 2 ) 2
24: δ ← max(δ (1) , δ (2) , δ (3) + δ (4) ) 25: return δ Algorithm 3 Validate a candidate solution of an integral equation Require: A polynomial integral operator K = (K ij ) 1 i,j p given by the (b ijk ) 1 i,j p 0 k κ ij , a polynomial right-hand side Ψ = (Ψ 1 , . . . , Ψ p ), a truncation order N val , (A, Λ) obtained from Algorithm 1 with Λ convergent to 0, and a candidate solution Φ • = (Φ • 1 , . . . , Φ • p ). Ensure: Two vectors of upper and lower bounds ε + and ε -such that Φ

• i - Φ i Ч 1 ∈ [ε - i , ε + i ] for 1 i p. All operations are performed in interval arithmetics 1: for k = 1 to p do P k ← Φ k + p j=1 K kj • Φ • j -Ψ k 2: for i = 1 to p do 3: Q i ← p k=1 A ik • п N val • P k + (1 -п n ) • P i 4: η i ← Q i Ч 1 5: end for 6: ε + ← (1 -Λ) -1 • η 7: for i = 1 to p do ε - i ← ((1 -D i • Λ) -1 • η) i 8: return ε + and ε - Estimating N val .
The following theorem (proved in Appendix A.3) provides a worst-case estimate for the minimal value of N val . Although theoretically interesting, this exponential bound is overpessimistic for a wide range of examples.

Theorem 4. Let B ij = κ ij k=0 b ijk Ч 1 and B = (B ij ) 1 i,j p . (i)
The following bound estimates the minimal possible value for N val making Algorithm 1 produce a contracting Newton-like operator:

N val = O dρ(B) 2 exp(2ρ(B)) ,
where ρ(B) denotes the spectral radius of B. (ii) For a given approximation Φ • of Φ and in order that Algorithm 3

computes error enclosures [ε - i , ε + i ] for ε i = Φ • i -Φ i Ч 1 with ε + i /ε - i
κ (for some κ > 1), N val for Algorithm 1 must be at least:

N val = O ν κ -1 d B 2 ∞ exp(2 B ∞ ) , with ν = max 1 i p ε ∞ /ε i , ε ∞ = max 1 i p ε i and B ∞ = max 1 i p p j=1 B ij the associated operator norm.

Example and Discussion

Consider the following order 1, two-dimensional system, for x ∈ [0, a] with a > 0, whose solutions (depicted in Figure 3) are highly oscillating functions.

Rescale it over [-1, 1] with the change of variable x = a 2 (1 + t): We give two different integral transforms associated to this equation. The integral transform described in [START_REF] Benoit | Rigorous uniform approximation of D-finite functions using Chebyshev expansions[END_REF] consists in integrating Equation [START_REF] Ortega | Iterative solution of nonlinear equations in several variables[END_REF] once, resulting into an integral equation for Y with polynomial kernel and righthand side given by: On the other side, the integral transform used in [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF] allows us to validate the derivative Φ = Y . The polynomial kernel and right-hand side are:

   y 1 = -x n y 2 y 2 = x m y 1 y 1 (0) = 1, y 2 (0) = 0 ⇒    Y 1 = -a 2 n+1 (1 + t) n Y 2 Y 2 = a 2 m+1 (1 + t) m Y 1 Y 1 (-1) = 1, Y 2 (-1) = 0 . ( 17
)
0 1 2 3 -1 0 1 x y 1 (x) y 2 (x)
K(t, s) = 0 a 2 n+1 (1 + s) n -a 2 m+1 (1 + s) m 0 , Ψ(t) = 1 0 . K(t,
K(t, s) = 0 a 2 n+1 (1 + t) n -a 2 m+1 (1 + t) m 0 , Ψ(t) = a 2 n+1 (1 + t) n 0 .
Now, K(t, s) is of degree 0 with respect to s, giving two polynomials b 001 and b 010 of respective degrees n and m.

Let's now focus on the first integral transform, with n = 5, m = 4, a = 3. Using the spectral method explained in Section 3.1 and implemented in our C library, we fix an approximation degree N app = 100 and obtain numerical approximations Y • 1 and Y • 2 , that must now be validated. The whole implemented procedure automatically computes and bounds for increasing values of N val the Newton-like operator T associated to the truncated operator K [N val ] . The approximate inverse is computed as an (2h , 2d ) almost-banded order 2(N val + 1) matrix. This process stops as soon as the total Lipschitz matrix returned by Algorithm 1 has a spectral radius less than 1. In case of failure of Algorithm 1, the procedure is relaunched with N val ← 2N val . For this example, we obtain N val = 1664, h = 48 and d = 304, giving the following Lipschitz matrix:

Λ = 9.73 • 10 -4 9.89 • 10 -2 3.60 • 10 -2 9.92 • 10 -2 , ρ(Λ) = 6.06 • 10 -2 .
The last step is performed by Algorithm 3. Given the numerical approximations

Y • 1 and Y • 2 , it computes η = Y • -T • Y • (Ч 1 ) 2
(the examples gives η 1 = 3.20 • 10 -3 and η 2 = 1.91 • 10 -3 ) and outputs the error enclosures given by Theorems 2 and 3:

ε - 1 = 2.99 • 10 -3 , ε + 1 = 3.41 • 10 -3 , ε - 2 = 1.78 • 10 -3 , ε + 2 = 2.04 • 10 -3 .
This whole process for this example takes about 30 seconds on a modern computer.

Comparison with decoupling/desingularization. In the case of polynomial coefficients, an alternative consists in decoupling the system to obtain p scalar LODEs of order p, at the cost of introducing singularities in the equations. As an example, the first component y 1 in [START_REF] Ortega | Iterative solution of nonlinear equations in several variables[END_REF] satisfies the following differential equation:

xy 1 -ny 1 + x n+m+1 y 1 = 0. ( 18
)
This equation is singular (its leading coefficient vanishes at 0), so our validation method cannot be used. However, with desingularization techniques [START_REF] Abramov | Apparent singularities of linear difference equations with polynomial coefficients[END_REF], one obtains a higher order but nonsingular equation, whose set of solutions (strictly) contains the ones of the singular equation. In our example, by differentiating Equation ( 18) n times and dividing the result by x:

y (n+2) 1 + 1 x d n dx n (x n+m+1 y 1 ) = 0. ( 19
)
By inverting the roles of n and m, one obtains a similar equation for y 2 . Hence, validating the approximation y of ( 17) can be done with the validation algorithm of [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF]. Several caveats must therefore be raised. Applying the integral operator of [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF] results into a totally intractable problem, since the minimal value for proving that T is contracting is far too large (in practice, we stopped at N val 10 6 ). This is due to the fact that this transform is used to validate the last derivative y (n+2) 1

, which increases very rapidly due to the highly oscillating behavior of y 1 . On the other hand, the integral transform of [START_REF] Benoit | Rigorous uniform approximation of D-finite functions using Chebyshev expansions[END_REF] yields a far more tractable problem: a truncation order N val = 750 is sufficient for our example. However, Equation ( 19) is very illconditioned because of the factorial terms created by the n differentiations. For instance, with classical double precision (53 bits), the scalar validation procedure is able to produce and bound a contracting Newton-like operator T (Algorithm 1), but Algorithm 3 outputs an upper bound ε + 1 = 2.57, which is 3 orders of magnitude larger than what was found with the vectorvalued validation method. The non D-finite case. In the case of nonpolynomial coefficients, there is no general method to decouple and desingularize the system. Moreover, these coefficients may not be known exactly, but only given as polynomial approximations together with rigorous error bounds. We believe that in such a general case, the vector-valued approach presented in this article is essential to approximate and validate the solution. For example, a successful application of our method to a station keeping problem of a satellite is given in [START_REF] Arantes Gilz | Validated Semi-Analytical Transition Matrix for Linearized Relative Spacecraft Dynam-ics via Chebyshev Polynomials[END_REF]. Future extensions include: validated expansions in other orthogonal polynomial bases for LODEs; automation and complexity analysis for some classes of nonlinear ODEs; formally proving this method in a proof assistant. Heartfelt acknowledgments are extended to Mioara Joldes , and Nicolas Brisebarre for their extraordinary support in the writing of this article, to Bruno Salvy for instructive discussions about D-finite functions, to Bogdan Pasca for his friendly culinary support.

This work was partially funded by the ANR FastRelax project.

Proof. Since A has nonnegative entries which bound those of B, it can be easily shown by the triangle inequality that for any exponent k 0, |B k ij | A k ij for all i, j ∈ 1, p . This directly implies the conclusion of Lemma 2. Proof of Lemma 1. First, 1-D i •Λ is nonsingular because D i •Λ is convergent to zero by use of Lemma 2, since its entries are clearly dominated by those of Λ ∈ M →0 p (R + ). Hence so is Λ -D i . Then we prove that 1 -Λ and 1 -D i • Λ both have positive determinant. The segment 1τ Λ (τ ∈ [0, 1]) connects 1 to 1 -Λ, and all these matrices are nonsingular, because τ Λ converges to zero according to Lemma 2. Since det(1) = 1 > 0, we get by connectedness that det(1 -Λ) > 0. A similar argument proves that det(1 -D i Λ) > 0, and hence det(D i -Λ) < 0.

Remember that for a nonsinglular matrix M , we have M -1 = (det M ) -1 Cof (M ) T , where Cof (M ) is the cofactor matrix of M , whose entries are the minors of M . Noticing that Cof (D i -Λ) ji = Cof (1 -Λ) ji for j ∈ 1, p and using the fact that det(D i -Λ) < 0, det(1 -Λ) > 0 and all entries in (1 -Λ) -1 are nonnegative, we conclude that all entries on the i-th row of (D i -Λ) -1 are non-positive.

A.2 Discussion About the Tightness of Error Enclosures of Section 2

In the one-dimensional case with a contracting operator of Lipschitz constant λ ∈ (0, 1), the ratio of the upper bound and the lower bound given by the Banach fixed-point theorem is equal to (1 + λ)/(1λ) > 1. This quantity does not depend on the approximation x • , and uniformly tends to 1 as λ → 0, justifying the principle: the more contracting the operator is, the tighter the obtained enclosure is. This section aims at extending this study to the vectorial case, for the bounds obtained from Theorems 2 and 3. The following lemma quantifies how much the obtained enclosure [ε - i , ε + i ] of ε i = d(x • , x ) i may overapproximate it.

Lemma 3. Let T be contracting of Lipschitz matrix Λ ∈ M →0 p (R + ), x its unique fixed-point, x • an approximation, ε = d(x • , x ) the approximation error, and η = d(x • , T • x • ). Then for all i ∈ 1, p :

ε + i := e T i • (1 -Λ) -1 • η e T i • (1 -Λ) -1 • (1 + Λ) • ε ε - i := e T i • (1 -D i • Λ) -1 • η e T i • (1 -D i • Λ) -1 • (1 + D i • Λ) • ε (21)
The overapproximation ratio can be bounded as follows:

ε + i ε - i 1 + αν i 1 -αν i , with ν i = ε ∞ ε i and α = 2 Λ ∞ 1 -Λ ∞ , (22)
provided that αν i < 1, where ε ∞ = max 1 i p |ε i | is the infinite norm over R p and Λ ∞ = max 1 i p p j=1 |Λ ij | the associated operator norm. In particular, Lemma 3 shows that the ratio now depends not only on Λ, but also on ε = d(x • , x ).

Proof. For the first inequality, we have:

ε + i = e T i • (1 -Λ) -1 • η e T i • (1 -Λ) -1 • (1 + Λ) • ε,
since η (1+Λ)•ε by Equation ( 5) and (1-Λ) -1 has nonnegative coefficients. For the second inequality, consider η = (1 + D i • Λ) • ε, so that:

η i = e T i • (1 -Λ) • ε η i , η j = e T j • (1 + Λ) • ε η j , for j = i
Since by Lemma 1, e T i • (1 -D i • Λ) -1 (the i-th line of (1 -D i • Λ) -1 ) has a positive entry on the i-th position and nonpositive entries on the other ones, we have:

ε - i = e T i • (1 -D i • Λ) -1 • η e T i • (1 -D i • Λ) • η.
In order to prove [START_REF] Van Den Berg | Rigorous numerics in dynamics[END_REF], we notice that:

e T i • (1 -Λ) -1 • (1 + Λ) • ε = ε i + 2e T i • n 1 Λ n • ε ε i + 2 n 1 Λ ∞ ε ∞ = ε i + α ε ∞ ,
under the condition that Λ ∞ < 1, which is automatic from assumption αν i < 1. The inequality:

e T i • (1 -D i • Λ) -1 • (1 + D i • Λ) • ε ε i -α ε ∞ , is proven similarly by considering D i • Λ instead of Λ and D i • Λ ∞ = Λ ∞ .
Taking the quotient of these two inequalities yields the bound [START_REF] Van Den Berg | Rigorous numerics in dynamics[END_REF].

Tightness cones. For a fixed Λ, the overapproximation ratio depends on the error distribution of the ε i . More specifically, using Inequality (21) of Lemma 3, the constraint ε + i /ε - i κ for some ratio κ > 1 is fulfilled when:

e T i • (1 -Λ) -1 • (1 + Λ) • ε κ • e T i (1 -D i • Λ) -1 • (1 + D i • Λ) • ε.
The intersection of these constraints for i ∈ 1, p defines a cone C κ in R p + of points ε for which the computed overapproximation ratio does not exceed converges to K, we may approximate:

(1 + K [N val ] ) -1 (Ч 1 ) p ≈ (1 + K) -1 (Ч 1 ) p = O (dB exp(2B)) .

• We therefore have: (ii) When a maximal overapproximation ratio κ > 1 is fixed, Lemma 4 provides a condition on the infinite norm Λ ∞ with Λ = DT (Ч 1 ) p :

Λ = O dB 2 exp(2B) N val , where Λ = (1 + K [N val ] ) -1 • (K -K [N val ] ) (Ч 1
Λ ∞ β 2 + β = O κ -1 ν with β = 1 ν κ -1 κ + 1
Again, we focus on the truncation error (since the approximation error can be made as small as desired for a given N val ) and we have Λ = (1 -

K [N val ] ) -1 • (K -K [N val ]
) (Ч 1 ) p . Similarly than for the spectral radius, we have:

Λ ∞ = O d B 2 ∞ exp(2 B ∞ ) N val .
We therefore obtain the desired estimate.

Figure 2 (

 2 b)). It is represented by the order p(N app + 1) square matrix depicted by blocks in Figure 2(c). By permuting the natural basis B p,Napp of (п Napp • Ч 1 ) p into B p,Napp : B p,Napp = (T 0 e 1 , . . . , T Napp e 1 , . . . . . . , T 0 e p , . . . , T Napp e p ), B p,Napp = (T 0 e 1 , . . . , T 0 e p , . . . . . . , T Napp e 1 , . . . , T Napp e p ),

Figure 2 :

 2 Figure 2: Almost-banded structure of integral operators

Figure 3 :

 3 Figure 3: Solution of (17) with n = 5, m = 4 and a = 3

  s), which is of degree 0 in t, is decomposed over the Chebyshev basis with respect to s into constant polynomials b 001 , b 101 , . . . , b n01 and b 010 , b 110 , . . . , b m10 .

Figure 4 :

 4 Figure 4: Tightness cones for the toy example

  ) p , and by taking the spectral radius (note that B and exp(2B) commute):ρ(Λ) = O dρ(B) 2 exp(2ρ(B)) N val ,which gives the estimate for N val to obtain a matrix with spectral radius less than 1.

http://perso.ens-lyon.fr/florent.brehard

Although commonly attributed to Perov[START_REF] Perov | On the Cauchy problem for a system of ordinary differential equations[END_REF] (in Russian), the idea of generalizing the Banach fixed-point theorem to generalized norms for investigating the componentwise errors in an iterative process first appeared in Kantorovich's work[START_REF] Kantorovich | Functional analysis in partially ordered spaces[END_REF] (in Russian).

Appendix

A.1 Complementary Proofs for Section 2 Proof of Theorem 2. (i) Endow X with the metric d ∞ (x, y) = d(x, y) ∞ = max 1 i p d i (x), so that (X, d ∞ ) is a complete metric space. For an order p square matrix A, define:

Then T k is a µcontraction for the d ∞ metric, so that the Banach theorem applies and gives x as the unique fixed-point of T k . Hence T can have at most one fixed point. From the following inequality:

Since Λ k → 0 as k → ∞, it is easy to prove that 1 -Λ is nonsingular, with nonnegative inverse (1 -Λ) -1 = k 0 Λ k 0. Therefore, multiplying both members of Inequality ( 20) by (1 -Λ) -1 is licit, so as to obtain the upper bound (4) for n = 0. The general bound for n 0 follows from the fact that T is Λ-Lipschitz.

Lemma 1. Let Λ ∈ M →0

p (R + ) be a convergent to zero nonnegative matrix. Then, for every i ∈ 1, p , Λ -D i is nonsingular and the entries on the i-th row of its inverse are nonnegative.

In order to prove the technical Lemma 1, we need the following fact:

p (R + ) a convergent to zero nonnegative matrix and B ∈ M p (R) a matrix whose entries are dominated by those of A:

for all i, j ∈ 1, p .

Then B is convergent to zero.

κ. Under a certain value for κ, C κ is empty, meaning that T is not contracting enough to achieve this ratio, whatever ε is. The cone C κ grows to a limit cone C ∞ as κ → +∞ : a point outside C ∞ means that the componentwise error distribution is so unbalanced that some lower bound ε - i is negative (hence rounded to zero). Figure 4 illustrates the cones C κ for different values of κ and the limit cone C ∞ arising in our toy example. Validation up to a given ratio. When the problem consists in finding an error enclosure [ε - i , ε + i ] for a given approximation x • of x such that ε + i /ε - i κ, the fixed-point validation operator T must be chosen sufficiently contracting. The following lemma characterizes how small the Lipschitz matrix Λ must be in that case.

then the upper and lower bounds computed using Theorems 2 and 3 satisfy

Proof. By Lemma 3, having:

. This in turn is equivalent to α β and finally Λ ∞ β/(2 + β).

A.3 Proof of Theorem 4

Proof. (i) The value of N val must be sufficiently large to ensure that the truncation error

) p is a convergent to zero matrix.

• We have as a direct consequence of the one-dimensional case [7, Lemma 3.5]:

• For i 0, the bound K i

is generalized from the one-dimensional case contained in the proof of [START_REF] Bréhard | Validated and numerically efficient Chebyshev spectral methods for linear ordinary differential equations[END_REF]Lemma 3.3], where C =


Hung Quoc

Hang Lu

Thi Nguyen

Nga Thuy

Tien Nguyen

Zung Nguyen

Tat Dat Tô
email: tatdat.to@gmail.com

Skin Lesion Analysis Towards Melanoma Detection for ISIC 2018

methods for the ISIC 2018 Competition: Skin Lesion Analysis Towards Melanoma Detection.

Introduction

In this paper we briefly discuss our methods for addressing Task 1 (lesion segmentation), Task 2 (lesion attribute detection) and Task 3 (classification) of the ISIC 2018 Competition on skin cancer.

For Task 1 and Task 2, we used a convolutional neuron network (CNN) model called Grass-Net, which we developed earlier this year for another problem (satellite image segmentation), and which is a combination of U-Net [START_REF] Ronneberger | Convolutional Networks for Biomedical Image Segmentation[END_REF] and ResNet [START_REF] He | Deep Residual Learning for Image Recognition[END_REF]. For Task 3, we use a new CNN model that we developed called Z-Net (Z for zig-zag), and pre-process the images using Task 1.

For Task 1 and Task 2, we used only the official data of the ISIC 2018 Competition [START_REF] Noel | Skin Lesion Analysis Toward Melanoma Detection: A Challenge at the[END_REF] and [START_REF] Tschandl | The HAM10000 Dataset: A Large Collection of Multi-Source Dermatoscopic Images of Common Pigmented Skin Lesions[END_REF]. For Task 3, due to the fact that the official dataset is very imbalanced (e.g., the number of images in the category "NV" is almost 60 times higher than "DF"), we added some data that we could find on the internet using searches by keywords, from free sources such as DermanetNZ, to reduce this imbalance.

2 Task 1: Lesion Segmentation

Model

We designed an architecture using Keras/Tensorflow, called GrassNet, which combines the U-Net of Ronneberger et al. [START_REF] Ronneberger | Convolutional Networks for Biomedical Image Segmentation[END_REF] and the ResNet structures of He et al. [START_REF] He | Deep Residual Learning for Image Recognition[END_REF]. More precisely, the network contains down-sampling and up-sampling layers as in the U-Net, and residual blocks are added to down-sampling layers as in ResNet. In addition, we concatenate down-sampling layers to up-sampling layers in order to preserve the information. Our GrassNet has 37 million trainable parameters.

The droupout with p = 0.5 is also added at the two last down-sampling convolution blocks to reduce overfitting. The Batch normalization and the Selu activation are added after each Convolution2D layer. The last activation function is Sigmoid.

The architecture operates on an input image of 256 × 256 pixels and produces a probability map of the same dimensions.

Preprocessing and Training

We resized each training image to 256 × 256, then trained the network using the Adam optimization algorithm [START_REF] Kingma | A Method for Stochastic Optimization[END_REF]. In order to improve the robustness of the model, we applied a set of transformations (augmentations): rotation with a random angle, flipping, and random noise.

Validation results

We obtained a Jaccard score of 0.781 on our first validation submission, after 280 epochs of training (11 minutes per epoch), on a PC with a GTX 1070 graphic card.

3 Task 2: Lesion Attribute Detection

Model

For Task 2, we use the same architecture GrassNet as for Task 1.

  Figure 1: GrassNet

Preprocessing

In order to improve accuracy, we first apply Task 1 to find the lesion on the image, cut down the image to a smallest possible square which contains the lesion, and then resize that square to 256 × 256, before feeding it to GrassNet.

Augmentation and Training

We trained each attribute separately, using the Adam optimization algorithm [START_REF] Kingma | A Method for Stochastic Optimization[END_REF]. Data augmentations similar to Task 1 (random rotations, flipping, and random noises) are used in order to improve the robustness of the model.

We did not submit our results for validation yet, due to lack of time. (We entered this competition late, on July 4th).

4 Task 3: Lesion Classification

Model

We designed a new network, called Z-Net (Z for zig-zag) by combining the GrassNet in Task 1 and an additional down-sampling structure. The idea is that the hidden segmentation inside Z-Net will improve the accuracy of the classification. Our test runs on this architecture show promising results compared to previously well-known architectures.

See Figure 2 for the schema of a version of Z-Net. A more sophisticated version also has bagging.

Preprocessing

We pre-process in a way similar to Task 2, using the results of Task 1.

Augmentation and Training

Since the data is extremely imbalanced, we use the up-sampling technique by multiplying the data with different ratios on different classes. Then we apply a random mix-up [START_REF] Zhang | mixup: Beyond Empirical Risk Minimization[END_REF] with random weights: for a random couple of image (img 1 , img 2 ) extracted from data set, we replace it by a new image img = α * img 1 + (1 -α) * img 2 with the same label as img 1 , where α is a random number in [0.5, 1]. Standard augmentation operations are also used (random rotations, random zooming, random noise, etc.) in the training process.