EXISTENCE OF TRAVELING WAVE SOLUTIONS FOR THE DIFFUSION POISSON COUPLED MODEL: A COMPUTER-ASSISTED PROOF

The Diffusion Poisson Coupled Model describes the evolution of a dense oxide layer appearing at the surface of carbon steel canisters in contact with a claystone formation. This model is a one dimensional free boundary problem involving drift-diffusion equations on the density of species (electrons, ferric cations and oxygen vacancies), coupled with a Poisson equation on the electrostatic potential and with moving boundary equations, which describe the evolution of the position of each unknown interfaces of the spatial domain. Numerical simulations suggest the existence of traveling wave solutions for this model. These solutions are defined by stationary profiles on a fixed size domain with interfaces moving both at the same velocity. In this paper, we present and apply a computer-assisted method in order to prove the existence of these traveling wave solutions. We also establish a precise and certified description of the solutions. Mathematics Subject Classification. 35C07, 35Q92, 47H10, 65G20, 65N35. Received December 18, 2020. Accepted July 13, 2021.

they show the efficiency of the developed methods and the relevance of the model. These numerical experiments have also highlighted the long-time behavior of the model: after a transient period, a kind of stationary regime is reached. It can be seen as a traveling wave solution: the size of the domain where the equations of the system are defined stays constant, both interfaces are moving with the same velocity while the densities of charge carriers and the electric potential have stationary profiles. The traveling wave solutions can be defined as solutions of a stationary DPCM, which will be also detailed in Section 2.
In [10], the existence of traveling waves solutions for the DPCM has been investigated. An existence result is established for a simplified model, where the electroneutrality in the oxide layer is assumed, so that the electric field is constant. Numerical methods for the computation of the traveling waves solutions are also proposed in [10] and numerical analysis of the simplified model is done.
The main novelty of the current paper is to prove the existence of traveling waves solution for the general DPCM, and to obtain a precise and certified description of the solutions (including the width of the oxide layer, and the value of the corrosion velocity). In order to do so, we use a computer-assisted argument, which allows us to validate a posteriori a numerically computed approximate solution. More precisely, we construct a fixed-point map which is based on a numerical solution, and use a combination of analytic estimates and interval arithmetic computations to prove that this map is contracting in an explicit neighborhood of the numerical solution. This yields both the existence of a true solution and guaranteed error bounds. For a broader overview of these computer-assisted arguments, we refer the reader to the surveys [13,14,18,22] and the books [16,21]. Our work follows in particular the techniques introduced in [15,19] for computer-assisted proofs using Chebyshev series.
The outline of the paper is the following. In Section 2, we present the evolutive DPCM and the associate stationary model which defines the traveling waves. The main result of the paper is given in this section and settled in Theorem 2.1. For the stationary model, we develop a spectral method for the computation of numerical solutions. This method is based on the expansion of the different unknowns into Chebyshev series. It leads to a nonlinear systems of equations which is solved by a Newton method. The numerical method is introduced in Section 3. Then, the aim is to certify the existence of a strong solution to the stationary DPCM in the neighborhood of a numerical solution. The tools needed for the proof are presented in Section 4. Finally, numerical experiments are given in Section 5. All the numerical results are computed by the spectral method and the existence of an exact solution in the neighborhood is certified. Appendix A is devoted to the presentation of the test case used in Section 5: values of the numerous physical parameters and definition of the associated scaled parameters.

The evolutive DPCM
The DPCM introduced in [4] describes the evolution of a dense oxide layer at the surface of carbon steel canisters in contact with a claystone formation under anaerobic conditions. As the size of the oxide layer is very thin compared to the waste overpack size, it is a one-dimensional model. The unknowns are the density of charge carriers -for the oxygen vacancies, for the electrons and for the ferric cations -, the electric potential and the position of the interfaces of the oxide layer: 0 and 1 . In this model 0 denotes the interface between the oxide layer and the claystone formation (outer interface) and 1 denotes the interface between the oxide layer and the carbon steel canisters (inner interface).
In this paper, we restrict our attention to the dimensionless model. The scaling process has been partly described in Section 5 of [9]. It will be detailed in Appendix A. When it is possible, we will use the generic notation for the carrier densities: can be either , or . For instance, we denote by the charge of each associated species: = 2, −1, 3 for = , , . We also denote by the current density and the mobility or diffusion coefficient for the corresponding species. As the time scale chosen for the scaling is the characteristic time scale of the ferric cations, the scaled quantities = / appear in the scaled system.
The equations for the carrier densities , , , as well as the boundary conditions, have the same form. For = , or , they are: − + ′ 0 ( ) = 0 ( ( 0 ( )), ( 0 ( ))), on = 0 ( ), ∀ > 0, (2.1b) − ′ 1 ( ) = 1 ( ( 1 ( )), ( 1 ( )), ), on = 1 ( ), ∀ > 0, (2.1c) The functions 0 and 1 are prescribed by the kinetics of the electrochemical reactions at the interfaces. They can be written in the following generic form: where the functions 0 , 1 , 0 and 1 are smooth and positive functions. We will specify their definitions in Appendix A. The electric potential satisfies the following Poisson equation: and the moving boundary equations are: Let us comment on the different parameters arising in the last equations: is the dimensionless applied potential between the metal and the claystone, see Figure 1 of [4]. hl is the net charge density of the ionic species in the host lattice.
-As the inner and outer interfacial structures behave like a capacitor, there are intrinsic voltage drops through these interfaces, called voltages of zero charge and denoted by Δ 0 and Δ 1 . -2 , 0 and 1 are positive dimensionless parameters coming from the scaling. -0 ( ) is the dissolution speed of the layer, given by 0 ( ) = 0 hl 0 ( 0 ( )) .
They all will be introduced in Appendix A.
In the system (2.1)-(2.3), is a given dimensionless applied potential. Let us already mention that it is deduced from a physical value (expressed in Volts and evaluated relatively to the electrode reference NHE) thanks to the scaling detailed in Appendix A. The system (2.1)-(2.3), where is given, is called the "potentiostatic case". In this case, one output quantity of interest is the total current of electrons at the inner interface defined by Then, it is also interesting to search = ( ) such that the electron charge balance at the inner interface is fixed to a given constant˜, which means tot =˜. (2.6) The system (2.1)-(2.3) with the additional unknown function of time and the additional equation (2.6) is referred as the "galvanostatic case". If˜= 0, we speak of free corrosion and is called "free corrosion potential". Finally, the system is supplemented with some initial conditions:

Reformulation of the evolutive DPCM on a fixed domain
In the one dimensional framework, it is always possible to define a change of variables that transform a system of partial differential equations written on a moving domain into a new system of partial differential equations defined on a fixed domain. It has been already done for the evolutive DPCM in [5], so that the numerical methods for the DPCM are defined on a fixed domain, with a fixed mesh. The same change of variables is used in [10] and in both cases the system is rewritten in [0, 1] × [0, ∞).
As we plan to expand the unknowns into Chebyshev series in order to use a spectral numerical method for the pseudo-stationary DPCM, we rewrite now the system of equations (2.1)-(2.3) in [−1, 1] × [0, ∞). Therefore, we use the following change of variable: It allows us to associate to every function ( = , , or ) defined on We also define the size of the domain by Therefore, we obtain We do not give here the details of the computation as they are classical and similar to those done in [5]. Moreover, we forget the bars and we come back to the notation instead of in the reformulated system. The equations on the carrier densities = , , then write: and are supplemented with the Fourier boundary conditions: The electric potential satisfies the following Poisson equation The moving boundary equations are written as: Finally, in the galvanostatic case, the relation (2.6) rewrites as

The stationary DPCM
Numerical experiments show the existence of traveling wave solutions for the DPCM (2.1)-(2.3), see [4]. These solutions do not depend on time and are defined in a fixed size domain, whose boundaries are moving at the same velocity. Therefore, they are solutions to the stationary form of the equations of the DPCM. Moreover, we set ( ) = ℓ the constant size of the domain and ′ 0 ( ) = ′ 1 ( ) = the constant velocity of the interfaces. The equations for the charge carrier densities = , , then write as: The equations for the electric potential write as: Finally, from the moving boundary equations (2.10), we deduce the following equations on the size ℓ and the velocity : Moreover, in the galvanostatic case, the additional equation (2.11) becomes We stress that the width ℓ of the oxide layer and the velocity of its interfaces are not input parameters, but part of the unknowns of the system, and that we have to solve for them. In this paper, we focus only on the potentiostatic case (2.12)-(2.14) to simplify the presentation, but the approach that we introduce could also handle the galvanostatic case. We both prove the existence of solutions of (2.12)-(2.14), and get quantitative, certified information about these solutions. These results will be presented in more details in Section 5.2, but Theorem 2.1 exhibits our main result for a given set of data. We emphasize that the parameter values pH = 8.5 and = 0.5 Volts play no particular role in our proof. We get similar results for different values of pH and in Section 5.2. There, we also give the corresponding values of and ℓ in physical units (remember that all the quantities in (2.12)-(2.14) were nondimensionalized).

Expansion into Chebyshev series and computation of a numerical solution
A numerical scheme, based on finite volumes, was introduced in [10] in order to obtain approximate solutions of the stationary DPCM. Here, we adopt a different strategy, which is more easily compatible with computerassisted proofs. The solutions (both approximate and rigorous) of (2.12)-(2.14) will be built as Chebyshev series. We start by recalling some basics facts about Chebyshev series, and we then introduce the numerical method we use for the computation of an approximate solution to (2.12)-(2.14). The extension to the galvanostatic case is straightforward: we simply need to incorporate the additional unknown and the additional equation (2.15), so we do not discuss it more in this paper.

Basics about Chebyshev series
We recall here a few results and notations about Chebyshev polynomials and series that are going to be useful in this work. For a more detailed exposition, further references and proofs, we refer to [20].

Functional analytic framework
In the sequel, we will need some functional analysis framework for the Chebyshev series. We will denote by any sequence ( ) ∈N in R N and throughout this note we will identify a function : [−1, 1] → R, at least Lipschitz continuous, with the sequence = ( ) ∈N of its Chebyshev coefficients.
For any sequence ∈ R N we define the -norm of as We also introduce the following Banach space Let us notice that for every , ∈ ℓ 1 we define * thanks to (3.6). This definition of a convolution product on ℓ 1 induces a natural Banach algebra structure on ℓ 1 , as stated in Lemma 3.5.
Proof. This follows from Definition 3.4 and the triangle inequality.

Reformulation of the stationary DPCM using Chebyshev series
The system (2.12)-(2.14) represents second order differential equations for the densities = , , and the electric potential . However, considering the currents for = , , and the electric field = as additional unknowns, it can be rewritten as a system of first order differential equations. More precisely, the inner equations (2.12a) and (2.13a) rewrite as: while the boundary conditions and the equations for the velocity and the thickness ℓ lead to: (3.10) Assume that the system (3.9) and (3.10) admits a smooth solution ( , , ( , ) = , , , , ℓ), so that , , and for = , , can be expanded into Chebyshev series: , for = , , . (3.11) We identify each function , , and with the sequence of its coefficients in the Chebyshev series: , , and . As the net charge density of the host lattice hl is a given constant, we can introduce the sequence hl ∈ R N defined by hl = ( hl , 0, . . .). Then, plugging the series expansions (3.11) into (3.9) and using the relation (3.6) and Proposition 3.1, we obtain the following infinite dimensional set of algebraic equations It is supplemented by the following relations, obtained by plugging the series expansions (3.11) into (3.10) and applying (3.1): Let us notice that since , = 0 for all ≥ 1 and = , , , the Chebyshev expansion of the function is given only by the first mode ,0 , i.e., = ( ,0 , 0, . . .). For better readability we will forget the subscript 0 and write . Moreover, in the sequel we will identify the real number with its natural injection in R N given by the sequence = ( , 0, . . .). We will also identify and ℓ with their natural injection in R N .

The formulation ( ) = 0
We now rewrite the infinite system of nonlinear equations (3.12) and (3.13) in the form ( ) = 0 where the unknown is and the function is defined by (3.14) where we recall that hl denotes the sequence given by hl = ( hl , 0, . . .).
We have rewritten the initial system of differential equations (2.12)-(2.14) as an infinite system of nonlinear equations of the form ( ) = 0. It remains now to prove that the existence of a solution to the system ( ) = 0 yields the existence of a smooth solution to the initial system (2.12)-(2.14). To this end we introduce an appropriate function space and in the sequel we will consider that the function given by (3.14) acts only on this space.
Let us first introduce I pot and I 1 pot the following index sets: and Definition 3.7. Let > 1 and ∈ (0, ∞) 10 be given. We define Remark 3.8. The norm || · || , depends on several parameters: ∈ (1, +∞) and ∈ (0, +∞) 10 , which must be carefully chosen in practice. We refer to Section 5.1 for a more in-depth discussion.

Computation of a numerical solution
From the infinite system of equations ( ) = 0, where is given by (3.14), we can easily deduce a finite nonlinear system of equations, just by truncating the Chebyshev modes of order ≥ 1 and higher, for a given . The unknowns of this new system are ( , , , , ) 0≤ ≤ −1 , , , , ℓ, and the size of the system is 5 + 5.
Let us denote by : ℓ 1 → R the finite dimensional projection obtained by truncating the Chebyshev modes of order We also denote by the natural injection from R 5 +5 to . We may now define We can compute an approximate solution, ∈ R 5 +5 to = 0 by solving numerically the finite dimensional problem [ ] = 0. We refer again to Section 5.1 for more details on the resolution of this finite dimensional problem. In the sequel, we use the same notation to denote ∈ R 5 +5 an approximate solution to [ ] = 0 and its injection into .

Presentation of the general strategy
We present in this section the strategy that will be used in order to prove the existence of a solution to (2.12)-(2.14). In Section 3, we have reformulated the problem as a zero finding problem ( ) = 0 for a suitable operator defined on the space . We are now going to introduce a Newton-like operator (see (4.6)) whose fixed points are in one-to-one correspondence with the zeros of . The existence and enclosure of the solution then follow by the contraction mapping theorem, once the operator is proven to be a contraction on some complete set. The following theorem, very reminiscent of the Newton-Kantorovich theorem, provides us with an efficient way of finding an explicit neighborhood of the numerical solution on which the operator is a contraction. Many similar versions of this theorem have been used in the last decades in computer-assisted proofs (see e.g. [1,11,16,27] and the references therein).

4)
where for any -linear ( ≥ 1) operator defined on we denote by ||| · ||| its operator norm. Define the radii polynomial as Assume that there exists > 0 such that ( ) < 0 and let and denote the two nonnegative roots of , with < . Then, provided < * , the operator : → defined as has a unique fixed point in ( , ) the closed ball of , centered at and of radius for all in [ min , max ) where min = and max = min Moreover since we assume that is an injective operator, then has a unique zero in ( , ) for all ∈ [ min , max ).
The proof simply consists in checking that is a contraction on ( , ) for all ∈ [ min , max ). We refer to the above-mentioned references for a detailed proof. -Since is simply a quadratic polynomial, the existence of an > 0 such that ( ) < 0 is equivalent to having and so these two conditions are sufficient conditions for to map ( , ) into itself. The restriction < + 2 is then enough to ensures that is contracting on this ball. All of this is only valid as long as ≤ * because of (4.4) (in practice, it is convenient to only have to control 2 in a small neighborhood of ).
-We are going to take for an approximate inverse of (︀ )︀ , and for † an approximation of (︀ )︀ itself (which will in fact be used to construct ). If we could take = i.e. a very strong contraction near . However, getting explicit estimates on (︀ (︀ )︀)︀ −1 can be very hard, which is why we introduce these approximations instead. The condition 0 + 1 < 1 tells us how good these approximations have to be. -Finally, if this condition 0 + 1 < 1 is satisfied, we only have to get good enough numerical approximation , or more precisely a small enough residual error , for the second condition 2 2 < (1 − ( 1 + 0 )) 2 to hold.
-In the sequel, we derive formulas for , 0 , 1 and 2 satisfying (4.1)-(4.4), which are explicit but cannot easily be evaluated by hand, since they depend on numerical data (and in particular on ). Therefore, we evaluate them with a computer, but using interval arithmetic (in our case with Intlab [17]) to ensure that the rounding errors are controlled.
Theorem 4.1 is the cornerstone of our computer-assisted proof. In Section 3.4, we have already introduced the function defined on the space ( , ‖ · ‖ , ) such that the solutions of = 0 correspond to the solutions of (2.12)-(2.14). Moreover in Section 3.5, we have defined ∈ R 5 × R 5 (identified with its injection in ) as an approximate solution of the finite dimensional problem [ ] = 0. Remark 4.3. Because and ℓ are not allowed to be equal to 0 in ( ( ) is not defined if or ℓ is equal to 0), ( , || · || , ) is not quite a Banach space. Nonetheless, any closed ball in which does not intersect the hyperplanes = 0 and ℓ = 0 is still a complete metric space, which is all we need to apply the contraction mapping theorem, so the conclusions of Theorem 4.1 are still valid, for all such that |¯| < and |l| < .
In order to use Newton-Kantorovich argument to prove Theorem 2.1, it remains: -to define the linear operators and † , -to define and compute the bounds , satisfying (4.1)-(4.4), -to check that ( ) given in (4.5) is negative for some > 0.
In what follows, we detail each of these steps.

Definition of the operators and †
Recalling that we want † to be an approximation of (︀ )︀ , we define for every ∈ the operator † as

Operator norms
In order to compute the bounds 0 , 1 and 2 in Theorem 4.1 we need to introduce some operator norms. First, let us consider a linear operator : ℓ 1 → ℓ 1 . We denote by ||| ||| the operator norm of , i.e.

||| ||| = sup
(4.10) It will be convenient to think of as an "infinite dimensional matrix", written in the canonical Schauder basis of ℓ 1 . That is, is characterized by the coefficients ( , ) , ∈N such that, for all ∈ ℓ 1 and all ∈ N, Similarly to the well know formula for matrix norms, we can express the operator norm of in terms of these coefficients. where ( , ) , ≥0 is the matrix representation of the operator .
In particular, we deduce that if the matrix has a finite number of non zero coefficients, then ||| ||| can be evaluated on a computer (and using interval arithmetic we can get a rigorous upper bound of this norm). It will be convenient to introduce the notation ·, = ( , ) ≥0 for every ≥ 0. Notice that we can then rewrite (4.11) as Now, let : → be a linear operator and let us consider the following block-representation of where the set I pot is given by (3.15). Due to the definition of , we note that for instance ( , ) : R → R, ( , ) : ℓ 1 → ℓ 1 for = , , and ( , ) : R → ℓ 1 . However, post-composing some blocks of by 1 , the natural injection from R to ℓ 1 , we can see each of these blocks as some linear operators from ℓ 1 or R to ℓ 1 . Moreover, this slight modification does not change the value of the operator norm of the blocks. For instance, ⃒ and we omit to write the composition by 1 in the sequel. Now, still considering : → , we slightly abuse the notation by applying the ℓ 1 operator norm component wise to , i.e. we define for , ∈ I pot , we use formula (4.12) to evaluate its operator norm.
Let us now introduce for ∈ (0, ∞) 10 the following weighted operator norm We recall that the practical choice of will be discussed in Section 5.1. If all the blocks of have a finite number of non zero coefficients, then the quantity ⃒ ⃒ ||| ||| ⃒ ⃒ can be evaluated on a computer (and rigorously upper-bounded using interval arithmetic). Moreover, we notice that Therefore, as soon as we can rigorously compute ⃒ ⃒ ||| ||| ⃒ ⃒ , we get a computable and rigorous upper bound for ||| ||| , .

The bound
We simply define as Since the vector has a finite number of non zero coefficients, i.e. = (0, 0, 0, 0, 0) for every ≥ , we have which implies that (︀ )︀ has a finite number of non zero coefficients. Moreover, since acts only diagonally on the tail of the elements of (︀ )︀ , (︀ )︀ also has a finite number of v non zero coefficients. Thus, the bound can be evaluated on a computer (using interval arithmetic).

The bound 0
Using the notations and definitions introduced in Section 4.3 we obtain the following result: satisfies We point out that, by construction of † and , each block ( , ) for , ∈ I pot has only a finite number of non zero coefficients (recall that the "tail" parts of † and act as the identity), and hence the ℓ 1 operator norms ⃒ can all be evaluated on a computer (using interval arithmetic).
Proof. The proof of the result follows directly from (4.13).

The bound 1
We now explain how to define a constant 1 such that First, we consider the linear operator = ( (︀ )︀ − † ) and we denote ( , ) for , ∈ I pot (recall definition (3.15) of I pot ) the block-representation of . Then, applying (4.13) we have We are going to bound the ℓ 1 operator norm of each ( , ) independently, using the following splitting: Then, we define for all , ∈ I pot where Γ in the case where the block ( , ) admits only one column. Let us notice, due to the structure of (︀ )︀ , † and , that Finally, we can define the constant 1 as which implies that 1 satisfies (4.17). The following result summarizes this approach and gives a precise definition of the constants Γ ( , ) tail for , ∈ I pot .
if the block ( , ) admits only one column. Moreover, let us introduce tail,2 are given in Table 1 for , ∈ I 1 pot and Γ we have where we recall that I 1 pot = { , , , , } ⊂ I pot and I 1 pot denotes the characteristic function of I 1 pot . Using the definition of (︀ )︀ and † we obtain ( = 0 for any 0 < < 2 − and any , ∈ I pot .
In particular, since we only consider here ≥ 2 , ( = 0 for any 0 < < , which yields ,0 ( Besides, by construction of , if ̸ = then = I 1 pot ( ) , as soon as ≥ or ≥ . Therefore, we get ,0 ( ,0 ( Moreover, still by construction of (︀ )︀ and † , we notice that ( for all ≥ 2 and ∈ I pot and ∈ I 1 pot , which implies Let us now give the explicit value of the first term in the right hand side in the case where = and = for instance. Then, we have Writing down the explicit value of the other terms, we observe, as in the former example, that this value only depends on the parity of ℓ. Thus, we obtain This yields definition (4.21) of Γ tail,2 given in Table 1. To this end we need to estimate the second term in the right hand side of the previous inequality. We consider first the case = and = . Then, we notice that Thus, remembering that we only consider here ≥ 2 , Using similar arguments we obtain the different values of Γ ( , ) tail,2 for = , , . Let us then consider the case = and = . We have , and a meticulous but rather straightforward analysis of the terms in ( ) leads to Then, we recall that = 0 for all ≥ , which allows us to rewrite the above sum as follows Applying similar arguments we obtain the values Γ ( , ) tail,2 and Γ ( , ) tail,2 given in Table 1. Finally, let us derive the bound Γ ( , ) tail,2 for = , , . To this aim, we notice that Then, a similar analysis that the one done to derive the bound Γ ( , ) tail,2 leads to where denotes the first element of the Schauder basis of ℓ 1 , i.e., = (1, 0, 0, . . .). Using ≥ 2 > , shifting the indices and applying the triangular inequality yields Recalling the definition (3.8) of the seminorm | · | , we get which concludes the proof of Proposition 4.6.

The bound 2
Proposition 4.7. Let > 1, ∈ (0, ∞) 10 , * > 0 and ∈ R 5 +5 (identified with its injection in ). Consider , † and defined by (3.14), (4.8) and (4.9) and the index set I pot given by (3.15) and define Then, Proof. We first notice that 2 ( ) is a bilinear operator. Similarly to what we have done in Section 4.3 for linear operator defined on , we use a block representation of bilinear operators defined on , this time with 3 indices:

Its operator norm, still denoted
and the following inequality holds Then we notice that for 1 , 1 , 2 ∈ I pot we have Applying the triangle inequality we obtain Finally, we conclude the proof of Proposition 4.7 using the fact that the -norm is a sub-multiplicative norm.
Let us notice that the computation of ⃒ for 1 , 2 ∈ I pot requires to take into account the tail part of ( 1, 2) for 1 = 2 ∈ { , , , , } (recall that by construction the tail part of the other blocks have all entries equal to zero). However, since it has a diagonal structure we can explicitly compute the operator norm of these blocks (using interval arithmetic). For instance, using ( , ) , = for ≥ and ∈ N, we have Regarding the second derivative of , most of the terms that appear are actually 0, or very straightforward to estimate, because the differential equations in the DPCM are merely quadratic in , , and . The only slightly more involved terms are the ones coming from the boundary conditions and the equation for the velocity (i.e. ( ) 0 , ( ) and ( ) in (3.14)), which are highly nonlinear, but the corresponding second derivatives can still be computed explicitly, and the supremum over the such that ⃦ ⃦ − ⃦ ⃦ , ≤ * is rigorously computed using interval arithmetic (in practice * is small).

Implementation details
The starting point of our theorem is an approximate solution to the stationary DPCM. As mentioned previously, we obtain such an approximate solution by applying Newton's method to the finite dimensional projection [ ] of .  Figure 2. We display here the evolution of the width ℓ of the oxide layer, of the corrosion speed and of the total current tot (all rescaled back to physical units) in terms of the potential (expressed in Volts and evaluated relatively to the electrode reference NHE), for several values of pH.
Of course, the initialization of Newton's method has to be done carefully. In order to get a suitable initial condition, we use a code from [10] which computes a solution of a simplified version of the stationary DPCM, where the coupling between the electric potential and the charge carriers is removed in the Poisson equation, i.e. with 0 in the r.h.s. of (2.13a). Starting from such an approximation, we use numerical continuation to gradually put the coupling back, until we get an approximate solution of (2.12)-(2.14). This is not the only option, and one could for instance integrate the time-dependent model for long enough, until we are close to the pseudo-stationary state, since it seems to be globally attracting.
The only thing that remains to be discussed before we present the results of the computer-assisted proofs is the important choice of the weights in the norm ‖ · ‖ , we use on the space . Actually, the crucial part for this problem is the careful choice of ∈ (0, ∞) 10 , while can be chosen a bit more carelessly. Indeed, must be strictly larger than 1 because we need terms like (4.21) to be small, and not too large because we do not want the various ‖ · ‖ norms appearing in the bounds to explode (this is related to the domain of analyticity of the functions , , and ), but we did not have to carefully select it in order for the proof to succeed: for all the results presented below = 1.1 is good enough. On the other hand, a naive choice for like = (1, 1, . . . , 1) never leads to a successful proof, because 1 ends up being way larger than 1 (remember (4.7)), so a deliberate finite (see Prop. 4.6), are independent of . Therefore, considering , 0 , 1 and 2 as function of , where the dependency in is explicit (see (4.14), (4.22) or (4.23)) it is cheap to numerically optimize for according to our needs. A possible optimization criteria, which has been successfully used in the past [6] and amounts to the computation of a Perron-Frobenius eigenvector, is to take such that 1 is minimal. However, for our current problem such a choice often leads to 2 being too large, and thus to the second condition in (4.7) no longer being satisfied. Therefore, we instead try to optimize for such that the two roots of (see (4.5)) are the furthest apart, under the constraint that 0 + 1 < 1. The optimization is done using an algorithm from Matlab's optimization toolbox. We emphasize that we do not actually care whether the we obtain is close to a global minimizer or not, as long as it is good enough for the conditions (4.7) to be satisfied. A slightly different approach to optimize the choice of the norm is discussed in [24].

Results
The set of parameters we use as test case is given in Appendix A. For three different values of pH, namely 7, 8.5 and 10, we numerically compute approximate solutions of (2.12)-(2.14) for 50 different values of the applied  potential regularly spaced between = 0 Volts and = 0.7 Volts. All these solutions are then rigorously validated using the procedure described in the paper, that is by evaluating the estimates derived in Section 4.4, and checking that the assumptions of Theorem 4.1 are satisfied. In Figure 2, we show the validated values of the width ℓ of the oxide layer, the velocity at which both interfaces move, and the total current tot (recall definition (2.11)), all scaled back to physical units. While we can (and did) compute and validate solutions up to = 0 Volts also when the pH is equal to 7 or 8.5, but the physical meaning of the obtained solutions is unclear, since the width of the oxide layer becomes smaller than 1 nanometer, and this is why the curves are truncated in these cases.
In Figures 3 and 4, we show some of the corresponding densities , , and (see also Fig. 1). Notice that, while the qualitative behavior of the profiles does not change much with , the interfaces get sharper when increases. We point out that our results match those obtained via the code CALIPSO [3].
The Matlab code used for this paper, including the implementation of all the bounds needed for the validation, can be found at [8].
Remark 5.1. While the setup presented in this paper only allows for the validation of solutions for a given set of parameters, we mention that a slight generalization of these techniques could be used to rigorously validate curves of solutions, when for instance all but one parameter is fixed (say all but ), and is varying, see e.g. [2,7,12,23,25,26]. 2) The DPCM system in physical variables has been introduced in [4]. The original densities of electrons, Fe 3+ cations and oxygen vacancies can be denoted , Fe and ox . The associate current densities will be denoted by , Fe and ox . The last unknowns are the electric potential Φ and the position of the interfaces 0 ( ) and 1 ( ).
the scaling on the densities and the potential implies that For the time, we use the scaling relative to the characteristic time of the cations. It means that the scale factor for the time is 2 0 / . This yields finally the convection-diffusion equations (2.1a) for the scaled densities = , , . Let us now focus on the boundary conditions for the densities in order to define the boundary functions ( 0 , 1 ) for = , , involved in (2.8). These boundary conditions are prescribed by the kinetics of the electrochemical reactions at the interfaces. At the interface oxide/solution, = 0 , the electrochemical reactions are the ferric release for the cations, the ferrous release and the proton reduction for the electrons and the oxygen exchange for the oxygen vacancies. At the interface oxide/metal, = 1 , they are the iron oxydation for the     cations, the electron exchange for the electrons and oxide host lattice growth for the oxygen vacancies. These boundary conditions can be written in the generic form, for = , Fe, ox, For the oxygen vacancies, we have: Fe .
Let us mention here that Fe 3+ is the activity of the ferric cations and Fe is the maximum occupancy for octahedral iron in cations in the oxide layer. For the electrons, the kinetics of the interface reactions yields: with Fe 2+ the activity of the ferrous cations, redox the redox potential in the solution, DOS the density of state of electrons in the metal.
Applying the scaling, we obtain that the functions 0 , 1 , 0 and 1 are defined by     , with the scaled kinetics coefficients defined in Table A.5.

A.3. Scaling of the moving boundary equations
Let us now finish with the parameters involved in the moving boundary equations (2.10). The Pilling-Bedworth ratio Π and the parameter depend on the molar volume of the metal Ω Fe . The dissolution kinetics contant 0 involved in the dissolution speed 0 is obtained after the scaling of the corresponding physical valuē 0 10 − pH . These parameters are presented in Table A.6.