Franck Davhys

Reval Langa

Morgan Pierre

Morgan Pierre Dedicated

Michel Pierre

A DOUBLY SPLITTING SCHEME FOR THE CAGINALP SYSTEM WITH SINGULAR POTENTIALS AND DYNAMIC BOUNDARY CONDITIONS

Keywords: Caginalp phase-field system, logarithmic potential, dynamic boundary conditions, maximal monotone operator, multivalued operator, convex splitting

research documents, whether they are published or not. The documents may come

Introduction

The Caginalp system has been proposed in [START_REF] Caginalp | An analysis of a phase field model of a free boundary[END_REF] to describe phase transition phenomena such as the melting-solidification in certain classes of materials. Dynamic boundary conditions have been introduced by physicists in the context of the Cahn-Hilliard equation to account for interactions of the bulk material with the walls [START_REF] Fischer | Novel surface modes in spinodal decomposition[END_REF][START_REF] Fischer | Diverging time and length scales of spinodal decomposition modes in thin flows[END_REF][START_REF] Fischer | Time-dependent density functional theory and the kinetics of lattice gas systems in contact with a wall[END_REF]. In such models, the time derivative appears explicitly in the boundary conditions.

The Caginalp system with dynamic boundary conditions is the initial and boundary value problem

               ∂ t w -∆w = -∂ t u, t > 0, x ∈ Ω, ∂ t u -∆u + f (u) -λu = w, t > 0, x ∈ Ω, ∂ t ψ -∆ Γ ψ + g(ψ) -αψ + ∂ n u = 0, t > 0, x ∈ Γ, ∂ n w| Γ = 0, u| Γ = ψ, w| t=0 = w 0 , u| t=0 = u 0 , ψ| t=0 = ψ 0 , (1.1)
where Ω denotes a bounded domain of R 2 or R 3 with sufficiently smooth boundary ∂Ω = Γ. In (1.1), λ and α are nonnegative constants, ∆ Γ is the Laplace-Beltrami operator and ∂ n is the outward normal derivative. The state variables (u, w) denote the order parameter and the temperature, respectively, and Ω is the domain which contains the material. A thermodynamically relevant choice for the nonlinearity f is the logarithmic function

f (s) = ln ( 1 + s 1 -s ) , s ∈ (-1, 1). (1.2)
The function g is assumed to be smooth and nonincreasing. More general assumptions on f and g with be given in (2.5)- (2.6).

The well-posedness and longtime behaviour of the Caginalp system with singular potentials and dynamic boundary conditions has been thoroughly analyzed in [START_REF] Cherfils | Corrigendum to: "Existence of global solutions to the Caginalp phase-field system with dynamic boundary conditions and singular potentials[END_REF][START_REF] Cherfils | Finite dimensional attractors for the Caginalp system with singular potentials and dynamic boundary conditions[END_REF][START_REF] Cherfils | Long time behavior of the Caginalp system with singular potentials and dynamic boundary conditions[END_REF][START_REF] Cherfils | On the Caginalp system with dynamic boundary conditions and singular potentials[END_REF] (see also [START_REF] Grasselli | The Caginalp phase-field system with coupled dynamic boundary conditions and singular potentials[END_REF]). The dynamic boundary condition for the order parameter is the third equation in (1.1). For the classical Caginalp system, this equation is usually replaced by the no-flux boundary condition for u. In the latter case, if f is the logarithmic potential (1.2), a separation property holds which guarantees that u stays away from the pure phases ±1 [START_REF] Cherfils | Some results on the asymptotic behavior of the Caginalp system with singular potentials[END_REF][START_REF] Grasselli | Long time behavior of solutions to the Caginalp system with singular potential[END_REF]. In contrast, for problem (1.2), the solution u may take the values ±1 on the boundary Γ, unless a sign condition holds for the nonlinearity on the boundary.

In the singular cases, solutions to (1.1) can no longer be interpreted in the usual sense. In order to remove the sign condition and to analyze the longtime behaviour of solutions, the authors used in [START_REF] Cherfils | Long time behavior of the Caginalp system with singular potentials and dynamic boundary conditions[END_REF] the notion of variational solution, by adapting to problem (1.1) the approach in [START_REF] Miranville | The Cahn-Hilliard equation with singular potentials and dynamic boundary conditions[END_REF] which was developed in the context of the Cahn-Hilliard equation. We refer the reader to the review papers [START_REF] Cherfils | The Cahn-Hilliard equation with logarithmic potentials[END_REF][START_REF] Miranville | The Cahn-Hilliard equation and some of its variants[END_REF] and references therein for more information on this subject.

From a numerical point of view, calculations for problem (1.1) with a logarithmic function were performed in [START_REF] Cherfils | Long time behavior of the Caginalp system with singular potentials and dynamic boundary conditions[END_REF] for regular solutions and for singular solutions, but only up to the singular time. An Allen-Cahn type problem which has some similarities with this problem (cf. Remark 4.5) was analyzed in [START_REF] Pierre | Global existence via a multivalued operator for an Allen-Cahn-Gurtin equation[END_REF], along with numerical computations of singular solutions in one space dimension. The numerical analysis and numerical computation of related problems involving dynamic boundary conditions and regular nonlinearities were also considered in, e.g., [START_REF] Antonietti | Discontinuous Galerkin approximation of linear parabolic problems with dynamic boundary conditions[END_REF][START_REF] Boyer | A DDFV method for a Cahn-Hilliard/Stokes phase field model with dynamic boundary conditions[END_REF][START_REF] Cherfils | A numerical analysis of the Cahn-Hilliard equation with nonpermeable walls[END_REF][START_REF] Cherfils | A numerical analysis of the Cahn-Hilliard equation with dynamic boundary conditions[END_REF][START_REF] Fukao | Structure-preserving finite difference schemes for the Cahn-Hilliard equation with dynamic boundary conditions in the one-dimensional case[END_REF][START_REF] Israel | Numerical analysis of a Cahn-Hilliard type equation with dynamic boundary conditions[END_REF][START_REF] Kovács | Numerical analysis of parabolic problems with dynamic boundary conditions[END_REF][START_REF] Nabet | Convergence of a finite-volume scheme for the Cahn-Hilliard equation with dynamic boundary conditions[END_REF]. Cahn-Hilliard type equations with logarithmic potential and classical boundary conditions have also drawn a lot of interest, cf., e.g. [START_REF] Bartels | Error control for the approximation of Allen-Cahn and Cahn-Hilliard equations with a logarithmic potential[END_REF][START_REF] Chen | Positivity-preserving, energy stable numerical schemes for the Cahn-Hilliard equation with logarithmic potential[END_REF][START_REF] Cherfils | Finite-dimensional attractors for the Bertozzi-Esedoglu-Gillette-Cahn-Hilliard equation in image inpainting[END_REF][START_REF] Cherfils | On the Bertozzi-Esedoglu-Gillette-Cahn-Hilliard equation with logarithmic nonlinear terms[END_REF][START_REF] Copetti | Numerical analysis of the Cahn-Hilliard equation with a logarithmic free energy[END_REF]. In these situations, one challenge is to adapt to the discretized problem the separation property which holds for the continuous time and space problem.

Up to now, the computation of singular solutions for Cahn-Hilliard type problems involving both a logarithmic potential and dynamic boundary conditions does not seem to have been addressed. Our purpose in this paper is to propose and to analyze a scheme which allows us to compute singular solutions to problem (1.1) even after the singularity occurs. A fundamental idea in our approach is to use the energy associated to the problem.

Indeed, let F and G denote an antiderivative of f and g, respectively. Then we have, assuming that (u, w) is a regular solution to (1.1),

d dt Ẽ(u(t), w(t)) + ∫ Ω |∇w(t)| 2 dx + ∫ Ω |∂ t u(t)| 2 dx + ∫ Γ |∂ t u(t)| 2 dσ = 0, (1.3)
where the energy Ẽ is defined by

Ẽ(u, w) def = ∫ Ω 1 2 |∇u| 2 + F (u) - λ 2 u 2 dx + ∫ Ω 1 2 |w| 2 dx + ∫ Γ 1 2 |∇ Γ u| 2 + G(u) - α 2 u 2 dσ.
Relation (1.3) shows in particular that Ẽ(u(t), w(t)) is nonincreasing as t increases. It can be obtained by multiplying the first equation in (1.1) by w, the second equation by ∂ t u, by summing the two resulting equations and integrating on Ω; an integration by parts using the dynamic boundary condition yields the result. Our time discretization is based on two ideas: a splitting in time, which decouples the resolution of the equations in the problem at each time step (as in, e.g. [START_REF] Batangouna | Convergence of exponential attractors for a time splitting approximation of the Caginalp phase-field system[END_REF]), and a convex splitting of the energy. Thus, the scheme is unconditionally uniquely solvable (Proposition 3.1). A discrete version of (1.3) is then obtained if the time step is small enough (Lemma 3.2).

By letting the time step tend to 0 and using a monotonicity argument, we show that the time semi-discrete solution converges to an energy solution of problem (1.1) (Theorem 3.4, our main result). This notion of solution is adapted from [START_REF] Grasselli | The Caginalp phase-field system with coupled dynamic boundary conditions and singular potentials[END_REF], where a problem similar to (1.1), with dynamic boundary conditions for the temperature as well, has been considered (our "energy solutions" are called "weak solutions" in [START_REF] Grasselli | The Caginalp phase-field system with coupled dynamic boundary conditions and singular potentials[END_REF]; they are more regular than the variational solutions described in [START_REF] Cherfils | Long time behavior of the Caginalp system with singular potentials and dynamic boundary conditions[END_REF]). Singularities are taken into account thanks to duality techniques involving a multivalued operator.

For the numerical computation of solutions, we use a finite element method for the space discretization of the problem. At each time step, a constrained convex minimization problem is solved in order to compute the order parameter u, with the constraint that u has values in [-1, 1]. For this purpose, the logarithmic nonlinearity f (cf. (1.2)) is regularized, thus making a gradient available for the energy.

Our paper is organized as follows. In Section 2, we introduce the functional framework and the notion of energy solution to problem (1.1). We describe the time semi-discrete scheme and its behaviour in Section 3. In Section 4, we focus on the analysis and numerical computation of 1d stationary singular solutions to (1.1). This particular situation was pointed out as a counter-example in [START_REF] Miranville | The Cahn-Hilliard equation with singular potentials and dynamic boundary conditions[END_REF]. Section 5 concludes the paper with numerical computations of regular and singular solutions to the Caginalp system in two space dimension.

Energy solutions

2.1. Main assumptions and notation. We set H := L 2 (Ω) and denote by (•, •) the scalar product in H (and also in H 2 and H 3 ) and by ∥ • ∥ the related norm. Next, we set V := H 1 (Ω) and denote by V ′ the (topological) dual of V . The duality between V ′ and V will be indicated by ⟨•, •⟩. Identifying H with H ′ through the scalar product of H, it is then well known that V ⊂ H ⊂ V ′ with continuous and dense inclusions. In other words, (V, H, V ′ ) is a Hilbert triplet (see, e.g., [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF]).

Since the system (1.1) also includes equations defined on Γ, we introduce some further spaces. Thus, we set H Γ := L 2 (Γ) and V Γ := H 1 (Γ) and denote by (•, •) Γ the scalar product in H Γ , by ∥ • ∥ Γ the corresponding norm, and by ⟨•, •⟩ Γ the duality between V ′ Γ and V Γ . In general, the symbol ∥ • ∥ X indicates the norm in the generic (real) Banach space X and ⟨•, •⟩ X stands for the duality between X and X ′ . We also denote by ∇ Γ the tangential gradient on Γ. We can thus define the spaces

H := H × H Γ and V := {z ∈ V : z| Γ ∈ V Γ } .
We introduce the H-scalar product in the following natural way, ( (u, ψ), (v, φ)

) H := (u, v) + (ψ, φ) Γ , and the associated norm is denoted ∥ • ∥ H . Next, we set, on V, ((u, v)) V := ∫ Ω (∇u • ∇v)dx + ∫ Γ u| Γ v| Γ + ∇ Γ u| Γ • ∇ Γ v| Γ dσ
and, analogously, we set, on V ,

((u, v)) V := ∫ Ω (∇u • ∇v)dx + ∫ Γ u| Γ v| Γ dσ.
Here and below, the restriction | Γ is understood in the sense of traces. We note that the norm ∥ • ∥ V related to the scalar product ((•, •)) V is equivalent to the usual one. It is not difficult to prove (see, e.g. [START_REF] Mugnolo | Dirichlet forms for general Wentzell boundary conditions, analytic semigroups, and cosine operator functions[END_REF]Lemma 2.1]) that the space V is dense in H. A characterization of the spaces V ′ and V ′ has been given in [START_REF] Grasselli | The Caginalp phase-field system with coupled dynamic boundary conditions and singular potentials[END_REF]Proposition 2.1].

We also define the continuous elliptic operators

A : V → V ′ , ⟨Av 1 , v 2 ⟩ := ∫ Ω ∇v 1 • ∇v 2 dx, A Γ : V Γ → V ′ Γ , ⟨A Γ ξ 1 , ξ 2 ⟩ Γ := ∫ Γ ∇ξ 1 • ∇ξ 2 dσ, A : V → V ′ , ⟨Av 1 , v 2 ⟩ V := ⟨Av 1 , v 2 ⟩ + ⟨A Γ ξ 1 , ξ 2 ⟩ Γ ,
where

ξ i = v i | Γ for i = 1, 2.
In particular, for v ∈ V , we have,

∥Av∥ V ′ def = sup ∥φ∥ V ≤1 ⟨Av, φ⟩ V = sup ∥φ∥ V ≤1 (∇v, ∇φ) ≤ ∥∇v∥, ( 2.1)
where we used the Cauchy-Schwarz inequality.

For v ∈ V ′ , we denote

⟨v⟩ := 1 |Ω| ⟨v, 1⟩,
where |Ω| is the measure of Ω. We notice that there exists a positive constant c Ω such that

|⟨v⟩| ≤ c Ω ∥v∥ V ′ , ∀v ∈ V ′ , (2.2) ∥v -⟨v⟩∥ V ′ ≤ c Ω ∥v∥ V ′ , ∀v ∈ V ′ , (2.3) ∥v -⟨v⟩∥ V ≤ c Ω ∥∇v∥, ∀v ∈ V.
(2.4)

Next, we state our hypotheses on the nonlinearities (cf. [START_REF] Cherfils | Long time behavior of the Caginalp system with singular potentials and dynamic boundary conditions[END_REF]). We assume that

         1. f ∈ C 2 (-1, 1), 2. f (0) = 0, lim s→±1 f (s) = ±∞, 3. f ′ (s) ≥ 0, lim s→±1 f ′ (s) = +∞, 4. f ′′ (s)sign(s) ≥ 0.
(2.5)

The logarithmic function f given by (1.2) satisfies these conditions. In (2.5), the condition f (0) = 0 is merely a normalization: in case f (0) ̸ = 0, then by changing f (s) into f (s) -f (0) and w into w -f (0) in system (1.1), we recover this normalization.

We assume that the (nonlinear) function g satisfies

g ∈ C 2 (R), g ′ (s) ≥ 0 and lim inf s→±∞ g ′ (s) ≥ κ 1 , (2.6)
for some κ 1 > 0. Note that we do not require g(0) = 0.

Remark 2.1. Let g be a function which belongs to C 2 ([-1, 1]). Then we can extend g to the whole real line to a C 2 function with compact support, also denoted g, and such that g(s) = g(s) -αs, s ∈ R, where g satisfies (2.6) and α > 0 is chosen large enough. Thus, we recover the assumptions on the nonlinearities f and g which were made in [START_REF] Cherfils | Long time behavior of the Caginalp system with singular potentials and dynamic boundary conditions[END_REF].

We set

F (s) = ∫ s 0 f (r)dr, s ∈ (-1, 1), G(s) = ∫ s 0 g(r)dr, s ∈ R.
If lim s→+1 F (s) exists, then we extend F by continuity to s = +1 and similarly, if lim s→-1 F (s) exists, we extend F by continuity to s = -1. This is the case, e.g. with the logarithmic function f (1.2). Moreover, we set F (s) = +∞ outside the effective domain of F ,

dom(F ) = {s ∈ R, F (s) < +∞} ⊂ [-1, 1].
Then, identifying f and g with maximal monotone graphs in R × R, we have f = ∂F and g = ∂G, ∂ representing the subdifferential of convex analysis (here in R) [START_REF] Attouch | Variational convergence for functions and operators[END_REF][START_REF] Barbu | Nonlinear semigroups and differential equations in Banach spaces[END_REF][START_REF] Brézis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF]. Next, we define the functional

J : H → R ∪ {+∞}, J(v, φ) := ∫ Ω F (v)dx + ∫ Γ G(φ)dσ,
where it is understood that one of, or both, the integrals in the expression of J may take the value +∞. We denote by ∂J the subdifferential of J in the space H, namely, for (u, ψ), (ξ, ζ) ∈ H, we have

(ξ, ζ) ∈ ∂J(u, ψ) def ⇐⇒ J(v, φ) ≥ J(u, ψ) + ( (ξ, ζ), (v, φ) -(u, ψ) ) H ∀(v, φ) ∈ H.
It is well known (cf, e.g., [3, Prop. 2.8 p. 67]) that the above condition is equivalent to ξ = f (u) a.e. in Ω and ζ = g(ψ) a.e. on Γ (2.7) (in particular, the function u then takes values in (-1, 1) a.e. in Ω).

We will also need a relaxation of the functional J. Thus, we introduce the restriction J V of J to the space V and its subdifferential ∂ V,V ′ J V in the V ′ -V duality. This, in view of the identification H ∼ H ′ , has to be understood as a maximal monotone operator in V × V ′ , namely, for U ∈ V and ξ ∈ V ′ , we have

ξ ∈ ∂ V,V ′ J V (U ) def ⇐⇒ J V ( Ũ ) ≥ J V (U ) + ⟨ξ, Ũ -U ⟩ V ∀ Ũ ∈ V.
(2.8)

A precise characterization of ∂ V,V ′ J in the spirit of (2.7) is not immediate to be given. In particular, it is reasonable to expect that ∂ V,V ′ J may be multivalued due to the bounded domain of f (see, e.g., [START_REF] Pierre | Global existence via a multivalued operator for an Allen-Cahn-Gurtin equation[END_REF]).

Existence and uniqueness of energy solutions.

Our assumptions on the initial data are

U 0 := (u 0 , ψ 0 ) ∈ V with J V (U 0 ) < +∞ and w 0 ∈ H. (2.9) Definition 2.2.
Let assumption (2.5)-(2.6) hold. Given (U 0 , w 0 ) which satisfies (2.9), we say that a pair (U, w) with U = (u, ψ) is an energy solution of problem (1.1) originating from (u 0 , ψ 0 , w 0 ) if

• U (0) = U 0 in H and w(0) = w 0 in H; • U ∈ C 0 ([0, +∞); H) ∩ L ∞ (0, T ; V) for any T > 0; • ∂ t U ∈ L 2 (0, T ; H) for any T > 0; • w ∈ C 0 ([0, +∞); H) ∩ L 2 (0, T ; V ) for any T > 0; • ∂ t w ∈ L 2 (0, T ; V ′ ) for any T > 0;
and the following relations hold for a.e. t in (0, +∞):

∂ t w(t) + ∂ t u(t) + Aw(t) = 0, in V ′ , (2.10) ∂ t U (t) + AU (t) + ξ(t) = (w(t), 0) + (λu(t), αψ(t)), in V ′ , (2.11) ξ(t) ∈ ∂ V,V ′ J V (U (t)), in V ′ .
(2.12)

We have:

Lemma 2.3. For two energy solutions (U i , w i ) departing from (U i 0 , w i 0 ), i = 1, 2, we have the following estimate on the difference ( Ũ , w) = (U 1 -U 2 , w 2 -w 2 ) in terms of the initial datum ( Ũ0 , w0 ) := (U 1 0 -U 2 0 , w 1 0 -w 2 0 ): ∥ Ũ (t)∥ 2 H + ∥ w(t)∥ 2 V ′ + ∫ t 0 ( ⟨A Ũ (s), Ũ (s)⟩ V + ∥ w(s)∥ 2
) ds

≤ Λe λ ′ t ( ∥ Ũ0 ∥ 2 H + ∥ w0 ∥ 2 V ′ ) , (2.13)
where Λ is a positive constant which is independent of t and of the initial data, and λ ′ = 2 max{λ, α} + 1.

Proof. We write (2.10) for (U 1 , w 1 ) and for (U 2 , w 2 ), we take the difference, and we integrate in time. We find that for all t ≥ 0,

w(t) + ũ(t) + A ∫ t 0 w(s)ds = w0 + ũ0 , in V ′ . ( 2.14)
Next, we write (2.11) for (U 1 , w 1 ) and for (U 2 , w 2 ), and we take the difference. This yields, for a.e. t ∈ (0, +∞),

∂ t Ũ (t) + A Ũ (t) + ξ 1 (t) -ξ 2 (t) = ( w(t), 0) + (λũ(t), α ψ(t)), in V ′ . (2.15)
We take the product of (2.14) by w(t) (in the V ′ -V duality) and the product of (2.15) by Ũ (t) (in the V ′ -V duality), we add the resulting equations, and we obtain 1 2

d dt ∫ t 0 ∇ w(s)ds 2 + 1 2 d dt ∥ Ũ (t)∥ 2 H + ∥ w(t)∥ 2 + ⟨A Ũ (t), Ũ (t)⟩ V +⟨ξ 1 (t) -ξ 2 (t), Ũ (t)⟩ V = λ∥ũ∥ 2 + α∥ ψ(t)∥ 2 Γ + ( w0 , w(t)) + (ũ 0 , w(t)
), for a.e. t ∈ (0, +∞). We used here that the term (ũ(t), w(t)) cancels. We write

( w0 , w(t)) = ( w0 , w(t) -⟨ w(t)⟩) + (⟨ w0 ⟩, w(t)) = 1 2 d dt ( 2 w0 , ∫ t 0 w(s) -⟨ w(s)⟩ds ) + (⟨ w0 ⟩, w(t)).
Using the monotonicity of ∂ V,V ′ J V and Young's inequality, we find

d dt Y(t) + ⟨A Ũ (t), Ũ (t)⟩ V + ∥ w(t)∥ 2 ≤ 2 max{λ, α}Y(t) + 2∥⟨ w0 ⟩∥ 2 + 2∥ũ 0 ∥ 2 , (2.16)
where

Y(t) = ∫ t 0 ∇ w(s)ds 2 - ( 2 w0 , ∫ t 0 w(s) -⟨ w(s)⟩ds ) + c∥ w0 ∥ 2 V ′ + ∥ Ũ (t)∥ 2 H .
Here, c > 0 is chosen large enough so that

Y(t) ≥ c 0 ( ∫ t 0 ∇ w(s)ds 2 + ∥ w0 ∥ 2 V ′ + ∥ Ũ (t)∥ 2 H ) (2.17)
for some c 0 > 0 (use (2.4) and Young's inequality). By (2.2), we have

∥⟨ w0 ⟩∥ 2 ≤ c ′ Ω ∥ w0 ∥ 2 V ′ . (2.18)
Using (2.18) in (2.16) and applying Gronwall's lemma to the resulting differential inequality, we obtain

Y(t) ≤ Ce λ ′ t ( ∥ Ũ0 ∥ 2 H + ∥ w0 ∥ 2 V ′ ) , t ≥ 0, (2.19)
with λ ′ = 2 max{λ, α} + 1. Using (2.14) and (2.1), we see that By using a regularization argument, we obtain the following result.

∥ w(t)∥ 2 V ′ ≤ C ′ ( ∫ t 0 ∇ w(s)ds 2 + ∥ w0 ∥ 2 V ′ + ∥ũ(t)∥ 2 + ∥ũ 0 ∥ 2 ) . ( 2
Theorem 2.5. Let (u 0 , ψ 0 , w 0 ) satisfy (2.9). Then the energy solution to problem (1.1) is also a variational solution (in the sense defined in [15, Definition 4.1]).

Proof. We first note [START_REF] Cherfils | Long time behavior of the Caginalp system with singular potentials and dynamic boundary conditions[END_REF]Theorem 4.3] that problem (1.1) possesses a unique variational solution for any initial condition which satisfies

(u 0 , ψ 0 , w 0 ) ∈ L ∞ (Ω) × L ∞ (Γ) × L 2 (Ω) and ∥u 0 ∥ L ∞ (Ω) ≤ 1, ∥ψ 0 ∥ L ∞ (Γ) ≤ 1.
This is a weaker requirement than (2.9). We approximate the initial datum (u 0 , ψ 0 , w 0 ) which satisfies (2.9) by a sequence (u k 0 , ψ k 0 , w k 0 ) such that for each k, (u k 0 , ψ k 0 , w k 0 ) is regular enough, satisfies (2.9), and

∥u k 0 -u 0 ∥ + ∥ψ k 0 -ψ 0 ∥ Γ + ∥w k 0 -w 0 ∥ → 0, as k → +∞.
We also introduce a C 1 regularization f N of the nonlinearity f , as in [START_REF] Cherfils | Long time behavior of the Caginalp system with singular potentials and dynamic boundary conditions[END_REF]. Since

f N is regular, the (unique) solution (u k N , ψ k N , w k N ) to (1.1)
where f is replaced by f N and where the initial condition is (u k 0 , ψ k 0 , w k 0 ), is regular. Now, we let first N tend to +∞, and then k tend to +∞. The proof of Theorem 4.3 in [START_REF] Cherfils | Long time behavior of the Caginalp system with singular potentials and dynamic boundary conditions[END_REF] shows that the limit (u, ψ, w) that we obtain is the variational solution to problem (1.1). On the other hand, the theory of maximal monotone operators [START_REF] Attouch | Variational convergence for functions and operators[END_REF][START_REF] Barbu | Nonlinear semigroups and differential equations in Banach spaces[END_REF][START_REF] Brézis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF] shows that for each k, (u k N , ψ k N , w k N ) tends to the energy solution (u k , ψ k , w k ) of (1.1) with initial condition (u k 0 , ψ k 0 , w k 0 ), as N tends to +∞. Lemma 2.3 implies that, as k tends to +∞, (u k , ψ k , w k ) tends to the energy solution of (1.1) which, by uniqueness of the limit, is (u, ψ, w).

From [15, Proposition 1], we deduce: Corollary 2.6. We assume that either lim s→±1 F (s) = +∞ or g(-1) + α < 0 < g(1) -α.

(2.21)

Then the energy solution satisfies -1 < u(x, t) < 1 for a.e. (x, t) ∈ Ω × R + and for a.e. (x, t) ∈ Γ × R + ; for all t 0 > 0 and for t ≥ t 0 , this solution solves (1.1) in the usual sense.

Condition (2.21) is known as the sign condition. For the logarithmic function f (1.2), we have lim s→±1 F (s) < +∞, and this sign condition is needed to guarantee that every solution to problem (1.1) is classical.

The time semi-discrete scheme

Let δt > 0 denote an arbitrary (fixed) time step. In a formal way (that is, assuming that the solution is regular enough), our time semi-discretization of (1.1) reads: let (u 0 , ψ 0 , w 0 ) = (u 0 , ψ 0 , w 0 ) be given and for n = 0, 1, 2 . . . , let (u n+1 , ψ n+1 , w n+1 ) solve

               (i) (u n+1 -u n )/δt -∆u n+1 + f (u n+1 ) = λu n + w n in Ω, (ii) (ψ n+1 -ψ n )/δt -∆ Γ ψ n+1 + g(ψ n+1 ) + ∂ n u n+1 = αψ n on Γ, (iii) u n+1 | Γ = ψ n+1 , (iv) (w n+1 -w n )/δt -∆w n+1 = -(u n+1 -u n )/δt in Ω, (v) ∂ n w n+1 | Γ = 0. (3.1)
This scheme is based on two types of splitting: a splitting in time, since we first solve the system (i)-(ii)-(iii) of (3.1), which gives us (u n+1 , ψ n+1 ), and secondly we solve (iv)-(v), which gives us w n+1 . We also use a convex splitting of the energy in (i)-(ii)-(iii), which ensures an unconditional unique solvability: the nonlinear (contractive) terms are treated implicitly, whereas the (expansive) terms λu and αψ are treated explicitly.

3.1. Unconditional unique solvability and energy estimate. The rigorous version of (3.1) reads as follows. Assume that (U 0 , w 0 ) = (U 0 , w 0 ) which satisfies (2.9) is given and for n = 0, 1, 2, . . . let (U n+1 , w n+1 ) ∈ V × V with U n+1 = (u n+1 , ψ n+1 ) be defined by

U n+1 minimizes Ū → E n ( Ū ) in V, (3.2) where E n ( Ū ) := 1 2δt ∥ Ū -U n ∥ 2 H + 1 2 ⟨A Ū , Ū ⟩ V + J V ( Ū ) -λ(u n , ū) -(w n , ū) -α(ψ n , ψ) Γ (3.3)
with Ū = (ū, ψ), and (once

U n+1 is computed) w n+1 ∈ V solves (w n+1 -w n )/δt + Aw n+1 = -(u n+1 -u n )/δt, in V ′ . (3.4)
We have: Proposition 3.1 (Unique solvability for all δt). Assume that (U 0 , w 0 ) = (U 0 , w 0 ) which satisfies (2.9) is given. Then there exists a unique sequence

(U n , w n ) n∈N in (V × H) N generated by the scheme (3.2)-(3.4). Moreover, for all n ∈ N, we have J V (U n ) < +∞.
Proof. Assume that (U n , w n ) satisfies U n ∈ V and w n ∈ H for some n ∈ N. Since J V is convex on V, the function E n , which is the sum of J V , of a continuous coercive quadratic form on V and of a continuous linear form on V, is strictly convex on V and lower semi-continuous. Moreover, by (2.5), we have F (s) ≥ 0 for all s ∈ dom(F ) and by (2.6), we also have

G(s) ≥ κ 1 4 s 2 -κ 2 , s ∈ R,
for some κ 2 ≥ 0. Since the linear terms in E n are dominated by the quadratic terms, we have

E n ( Ū ) ≥ c 0 ∥ Ū ∥ 2 V -c 1 , Ū ∈ V,
for some c 0 > 0 small enough and some c 1 ≥ 0. Thus, E n ( Ū ) → +∞ as ∥ Ū ∥ V → +∞ and so E n has a unique minimizer U n+1 in the Hilbert space V (see, e.g., [START_REF] Ekeland | Convex analysis and variational problems[END_REF]). Moreover, E n (U n+1 ) ≤ E n (U n ) < +∞, so J V (U n+1 ) < +∞ as well. Finally, the Lax-Milgram theorem shows that problem (3.4) has a unique solution w n+1 ∈ V .

We define the energy associated to U through

E( Ū ) := 1 2 ⟨A Ū , Ū ⟩ V + J V ( Ū ) - λ 2 ∥ū∥ 2 - α 2 ∥ ψ∥ 2 Γ , Ū = (ū, ψ) ∈ V. (3.5)
Lemma 3.2 (Energy estimate for small δt). Assume that δt ≤ 1/2 and let

(U n , w n ) ∈ V × H such that J V (U n ) < +∞. Then the solution (U n+1 , w n+1 ) of (3.2)-(3.4) satisfies E(U n+1 ) + 1 2 ∥w n+1 ∥ 2 + 1 4δt ∥U n+1 -U n ∥ 2 H + δt∥∇w n+1 ∥ 2 ≤ E(U n ) + 1 2 ∥w n ∥ 2 . (3.6) Proof. By (3.2), we have E n (U n+1 ) ≤ E n (U n ), which reads 1 2δt ∥U n+1 -U n ∥ 2 H + 1 2 ⟨AU n+1 , U n+1 ⟩ V + J V (U n+1 ) +λ(u n , u n -u n+1 ) + α(ψ n , ψ n -ψ n+1 ) Γ -(w n , u n+1 -u n ) ≤ 1 2 ⟨AU n , U n ⟩ V + J V (U n ),
where U n = (u n , ψ n ). Using the well-known identity

(a, a -b) = 1 2 ∥a∥ 2 - 1 2 ∥b∥ 2 + 1 2 ∥a -b∥ 2 (3.7)
(and similarly for the scalar product in H Γ ) in the terms involving λ and α, we find

1 2δt ∥U n+1 -U n ∥ 2 H + 1 2 ⟨AU n+1 , U n+1 ⟩ V + J V (U n+1 ) -(w n , u n+1 -u n ) - λ 2 ∥u n+1 ∥ 2 + λ 2 ∥u n+1 -u n ∥ 2 - α 2 ∥ψ n+1 ∥ 2 Γ + α 2 ∥ψ n+1 -ψ n ∥ 2 Γ ≤ 1 2 ⟨AU n , U n ⟩ V + J V (U n ) - λ 2 ∥u n ∥ 2 - α 2 ∥ψ n ∥ 2 Γ . (3.8)
Next, we write

-(w n , u n+1 -u n ) = -(w n+1 , u n+1 -u n ) + (w n+1 -w n , u n+1 -u n ) = (w n+1 , w n+1 -w n ) + δt∥∇w n+1 ∥ 2 + (w n+1 -w n , u n+1 -u n ) = 1 2 ∥w n+1 ∥ 2 - 1 2 ∥w n ∥ 2 + 1 2 ∥w n+1 -w n ∥ 2 +δt∥∇w n+1 ∥ 2 + (w n+1 -w n , u n+1 -u n )
where, in the second line, we used (3.4) and in the third line, we used (3.7). By the Cauchy-Schwarz inequality and Young's inequality, we have

|(w n+1 -w n , u n+1 -u n )| ≤ 1 4δt ∥u n+1 -u n ∥ 2 + δt∥w n+1 -w n ∥ 2 ,
and so (3.8) implies

1 4δt ∥U n+1 -U n ∥ 2 H + E(U n+1 ) + 1 2 ∥w n+1 ∥ 2 + δt∥∇w n+1 ∥ 2 + ( 1 2 -δt ) ∥w n+1 -w n ∥ 2 ≤ E(U n ) + 1 2 ∥w n ∥ 2 .
(3.9)

Since δt ≤ 1/2, this yields (3.6).

Remark 3.3. The energy estimate (3.6) is a discrete version of (1.3).

3.2. Convergence as δt → 0. For a time step δt > 0, let (U n , w n ) n≥0 be a sequence generated by the time semi-discrete scheme (3.2)-(3.4). We define the following functions from R + into V:

U δt (t) = ((n + 1) -t/δt)U n + (t/δt -n)U n+1 , t ∈ [nδt, (n + 1)δt) (n ∈ N), U δt (t) = U n+1 , t ∈ [nδt, (n + 1)δt) (n ∈ N), U δt (t) = U n , t ∈ [nδt, (n + 1)δt) (n ∈ N).
We define similarly the functions w δt , w δt and w δt from R + to H associated to the sequence (w n ) n≥0 . Note that

∂ t U δt = U n+1 -U n δt in D ′ ( (0, +∞); V ) (3.10)
and

∂ t w δt = w n+1 -w n δt in D ′ ( (0, +∞); H ) . (3.11)
Here, as usual, the notation D ′ means that the equality holds in the sense of distributions (see, e.g., [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF]).

Theorem 3.4. Let (U 0 , w 0 ) = (U 0 , w 0 ) such that (2.9) holds. Then the solution (U δt , w δt ) associated to scheme (3.2)-(3.4) tends to the energy solution (U, w) of problem (1.1) in the following sense, as δt → 0:

U δt → U weakly-⋆ in L ∞ (0, +∞; V), U δt → U strongly in C 0 ([0, T ]; H), for any T > 0, ∂ t U δt → ∂ t U weakly in L 2 (0, +∞; H), w δt → w weakly-⋆ in L ∞ (0, +∞; H), w δt → w strongly in C 0 ([0, T ], for any T > 0, ∂ t w δt → ∂ t w weakly in L 2 (0, T ; V ′ ), for any T > 0.
Proof. By Proposition 3.1, for all n we have J V (U n ) < +∞, so that -1 ≤ u n ≤ 1 for a.e. x ∈ Ω and -1 ≤ ψ n ≤ 1 for a.e. x ∈ Γ.

(3.12)

By induction, estimate (3.6) yields

E(U n ) + 1 2 ∥w n ∥ 2 + 1 4δt n-1 ∑ k=0 ∥U k+1 -U k ∥ 2 H + δt n-1 ∑ k=0 ∥∇w k+1 ∥ 2 ≤ E(U 0 ) + 1 2 ∥w 0 ∥ 2 ,
for all n ≥ 1. Using (3.12) and the definition (3.5) of E, this yields

1 2 ⟨AU n , U n ⟩ + J V (U n ) + 1 2 ∥w n ∥ 2 + 1 4δt n-1 ∑ k=0 ∥U k+1 -U k ∥ 2 H +δt n-1 ∑ k=0 ∥∇w k+1 ∥ 2 ≤ E(U 0 ) + 1 2 ∥w 0 ∥ 2 + λ 2 |Ω| 1/2 + α 2 |Γ| 1/2 , ( 3.13)
for all n ≥ 1. Since the right-hand side above is finite by assumption, this shows that (U n ) n is bounded in V and (w n ) n is bounded in H. Thus, (U δt ), (U δt ), (U δt ) are bounded in L ∞ (0, +∞; V) and (w δt ), (w δt ), (w δt ) are bounded in L ∞ (0, +∞; H). By letting n tend to +∞ in (3.13), we see that

1 4δt +∞ ∑ k=0 ∥U k+1 -U k ∥ 2 H + δt +∞ ∑ k=0 ∥∇w k+1 ∥ 2 ≤ C, ( 3.14)
where C < +∞ is the right-hand side of (3.13). Thanks to (3.10), this can be rewritten 1 4

∫ +∞ 0 ∥∂ t U δt ∥ 2 H dt + ∫ +∞ 0 ∥∇w δt ∥ 2 dt ≤ C, (3.15) which shows that (∂ t U δt ) is bounded in L 2 (0, +∞; H) and (∇w δt ) is bounded in L 2 (0, +∞; H d ) (with d = 2 if Ω ⊂ R 2 and d = 3 if Ω ⊂ R 3 ). Moreover, δt +∞ ∑ k=0 ∥U k+1 -U k ∥ 2 H = ∫ +∞ 0 ∥U δt -U δt ∥ 2 H dt ≤ 4δt 2 C, (3.16)
and similarly,

∫ +∞ 0 ∥U δt -U δt ∥ 2 H dt ≤ 4δt 2 C (3.17)
From the estimates above, we deduce that, up to a subsequence, (U δt ), (U δt ) and (U δt ) converge weakly-⋆ in L ∞ (0, +∞; V) to the same limit U , as δt tends to 0. The limit U is the same in L 2 (0, +∞; H) by (3.16)-(3.17), and therefore also in D ′ (0, +∞; H). Moreover, ∂ t U ∈ L 2 (0, +∞; H) by (3.15). By a well-known result [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF], for any T > 0, the space

{ Ũ ∈ L ∞ (0, T ; V), ∂ t Ũ ∈ L 2 (0, T ; H) } is compactly imbed-
ded in the space C 0 ([0, T ]; H). Thus, (U δt ) converges to U strongly in C 0 ([0, T ]; H), for any T > 0. In particular, the initial condition U (0) = U 0 is satisfied. Using (3.4) and (2.1), we have

∥w n+1 -w n ∥ V ′ ≤ δt∥∇w n+1 ∥ + ∥u n+1 -u n ∥ H . Thus, δt +∞ ∑ n=0 ∥w n+1 -w n ∥ 2 V ′ ≤ 2δt 3 +∞ ∑ n=0 ∥∇w n+1 ∥ 2 + 2δt +∞ ∑ n=0 ∥u n+1 -u n ∥ 2 H . (3.18)
Using the estimates above, we see that

∫ +∞ 0 ∥w δt -w δt ∥ 2 V ′ dt ≤ ∫ +∞ 0 ∥w δt -w δt ∥ 2 V ′ dt ≤ 8Cδt 2 .
From these estimates, we deduce that, up to a subsequence, (w δt ), (w δt ), and (w δt ) converge weakly-⋆ in L ∞ (0, +∞; H) and weakly in L 2 (0, T ; V ) for any T > 0 to the same limit w.

It remains to prove that the limit, which has been denoted (U, w) for notational convenience, is indeed the energy solution to problem (1.1).

First, owing to (3.10)-(3.11), equation (3.4) can be written

∂ t w δt = -Aw δt -∂ t u δt in L 2 (0, T ; V ′ ), for any T > 0. (3.19)
Thus, (∂ t w δt ) is bounded in L 2 (0, T ; V ′ ) for any T > 0 and letting δt tend to 0 in this (linear) equation, we find that

∂ t w = -Aw -∂ t u in L 2 (0, T ; V ′ ), for any T > 0. (3.20)
Since, for any T > 0, the space [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF], we also obtain that (w δt ) converges strongly to w in C 0 ([0, T ]; H), for any T > 0. In particular, the initial condition w(0) = w 0 is satisfied.

{ w ∈ L 2 (0, T ; V ), ∂ t w ∈ L 2 (0, T ; V ′ ) } is compactly imbedded in the space C 0 ([0, T ]; H)
Secondly, owing to the definition (2.8) of ∂ V,V ′ J V , the Euler-Lagrange equation of the minimization problem (3.2)-(3.3) reads

1 δt (U n+1 -U n ) + AU n+1 + ξ n+1 = (λu n , αψ n ) + (w n , 0) in V ′ , with ξ n+1 ∈ ∂ V,V ′ J V (U n+1
). This can be rewritten

∂ t U δt + AU δt + ξ δt = (λu δt , αψ δt ) + (w δt , 0) in L 2 (0, T ; V ′ ), (3.21)
for any T > 0, where

ξ δt (t) def = ξ n+1 ∈ ∂ V,V ′ J V (U δt (t)) for all t ∈ [nδt, (n + 1)δt) (n ∈ N). (3.22)
We set T > 0 and we introduce the maximal monotone operator

B T in L 2 (0, T ; V)× L 2 (0, T ; V ′ ) defined for all (U, ξ) ∈ L 2 (0, T ; V) × L 2 (0, T ; V ′ ) by [7, Exemple 2.3.3. p.25 ] ξ ∈ B T (U ) ⇐⇒ ξ(t) ∈ B T U (t) for a.e. t ∈ (0, T ).
Since B T is a maximal monotone operator [7, p.22], for all (U, ξ) ∈ L 2 (0, T ; V) × L 2 (0, T ; V ′ ), we have

ξ ∈ B T (U ) ⇐⇒ ∫ T 0 ⟨ξ(t) -ξ(t), U (t) -Ũ (t)⟩ V dt ≥ 0 (3.23) for all ( Ũ , ξ) ∈ L 2 (0, T ; V) × L 2 (0, T ; V ′ ) such that ξ ∈ B T ( Ũ ). Let ( Ũ , ξ) ∈ L 2 (0, T ; V) × L 2 (0, T ; V ′ ) such that ξ ∈ B T ( Ũ ). From (3.22)-(3.

23), we deduce that

I δt def = ∫ T 0 ⟨ξ δt -ξ, U δt -Ũ ⟩ V dt ≥ 0. ( 3
I δt = ∫ T 0 ⟨-∂ t U δt -AU δt + (λu δt , αψ δt ) + (w δt , 0) -ξ, U δt -Ũ ⟩ V dt.
We write I δt = I 1 δt + I 2 δt + I 3 δt where

I 1 δt = - ∫ T 0 ⟨∂ t U δt , U δt -Ũ ⟩ V dt, I 2 δt = - ∫ T 0 ⟨AU δt , U δt -Ũ ⟩ V dt, I 3 δt = ∫ T 0 ⟨(λu δt , αψ δt ) + (w δt , 0) -ξ, U δt -Ũ ⟩ V dt.
Arguing as in (3.17), we see that

∫ T 0 ∥U δt -U δt ∥ 2 H dt ≤ 4δt 2 C.
Using that (U δt ) converges strongly to U in C 0 ([0, T ]; H), we obtain that (U δt ) converges strongly to U in L 2 (0, T ; H). Thus, ((λu δt , αψ δt )) converges strongly to (λu, αψ) in L 2 (0, T ; V ′ ). Similarly, (w δt ) converges strongly to w in L 2 (0, T ; H). This shows that

I 3 δt → ∫ T 0 ⟨(λu, αψ) + (w, 0) -ξ, U -Ũ ⟩ V dt, (3.25)
as δt → 0. Next, we write

I 1 δt = - ∫ T 0 ⟨∂ t U δt , U δt ⟩ V dt - ∫ T 0 ⟨∂ t U δt , U δt -U δt ⟩ V dt + ∫ T 0 ⟨∂ t U δt , Ũ ⟩ V dt.
In the right-hand side above, the last integral tends to ∫ T 0 ⟨∂ t U, Ũ ⟩ V dt and the second integral tends to 0 thanks to the Cauchy-Schwarz inequality and the estimates (3.15) and (3.17). The first integral reads (see, e.g. [39, Lemma 1.2, Ch. III])

∫ T 0 ⟨∂ t U δt , U δt ⟩ V dt = 1 2 ∥U δt (T )∥ 2 H - 1 2 ∥U δt (0)∥ 2 H ,
and this tends to

1 2 ∥U (T )∥ 2 H - 1 2 ∥U (0)∥ 2 H = ∫ T 0 ⟨∂ t U, U ⟩ V dt,
owing to the strong convergence of (U δt ) in C 0 ([0, T ], H). Thus,

I 1 δt → - ∫ T 0 ⟨∂ t U, U -Ũ ⟩ V dt, (3.26)
as δt → 0. The term I 2 δt reads

I 2 δt = - ∫ T 0 ⟨AU δt , U δt ⟩ V dt + ∫ T 0 ⟨AU δt , U ⟩ V dt.
The second integral in the right-hand side above tends to ∫ T 0 ⟨AU, U ⟩ V dt, and by lower semi-continuity of the semi-norm in the Hilbert space L 2 (0, T ; V), we have

∫ T 0 ⟨AU, U ⟩ V dt ≤ lim inf δt→0 ∫ T 0 ⟨AU δt , U δt ⟩ V dt. Thus, - ∫ T 0 ⟨AU, U -Ũ ⟩ V dt ≥ lim sup δt→0 I 2 δt . (3.27)
Summing up (3.24)-(3.27), we have proved that

∫ T 0 ⟨-∂ t U -AU + (λu, αψ) + (w, 0) -ξ, U -Ũ ⟩ V dt ≥ 0.
From (3.23), we deduce that ξ ∈ B T (U ) where

ξ = -∂ t U -AU + (λu, αψ) + (w, 0) in L 2 (0, T ; V ′ ).
Recalling the definition of B T , this reads

∂ t U (t) + AU (t) + ξ(t) = (λu(t), αψ(t)) + (w(t), 0) in V ′ , with ξ(t) ∈ ∂ V,V ′ J(U (t)
), for a.e. t ∈ (0, T ). This relation is valid for any T > 0. It shows that (U, w) is the energy solution to problem (1.1). By uniqueness of the energy solution, the whole sequence (U δt , w δt ) tends to (U, w). The proof is complete.

Analysis and numerical approach for a singular 1d stationary case

In this section, we focus on the following boundary value problem in (-L, L) (0 < L < +∞): for a given K ≥ 0, find the function y which solves

{ -y ′′ + f (y) = 0 in (-L, L), y ′ (-L) = K, y ′ (+L) = K, ( 4.1)
where f is the (singular) logarithmic potential

f (s) = ln ( 1 + s 1 -s ) , s ∈ (-1, 1). (4.2)
This problem was pointed out as a counter-example in [START_REF] Miranville | The Cahn-Hilliard equation with singular potentials and dynamic boundary conditions[END_REF]. It can be seen as a stationary case of the Caginalp system (1.1) for Ω = (-L, L) by choosing w = 0, y = u, λ = α = 0 and g(s) = -K if x = +L, g(s) = +K if x = -L (since Ω is an interval, the Laplace-Beltrami operator does not appear in this case).

4.1. Analysis. For small K, problem (4.1) has a (unique) classical solution, but for large K, this is no longer the case (cf. Proposition 4.2) and it is necessary to introduce another notion of solution (note that for K ≤ 0, the situation is similar since -y is a solution of the ODE). In the spirit of our scheme, we consider energy solutions.

More precisely, assume that y is a classical solution of (4.1). Then, on multiplying the first equation by a test function φ and integrating by parts, we find that y satisfies

∫ L -L y ′ (x)φ ′ (x)dx -Kφ(+L) + Kφ(-L) + ∫ L -L f (y(x))φ(x)dx = 0, (4.3)
for all φ ∈ C 1 ([-1, 1]). Thus, y is a critical point of the functional

E K (z) := ∫ L -L (z ′ ) 2 2 + F (z)dx -Kz(+L) + Kz(-L), (4.4)
where

F (s) = ∫ s 0 f (σ)dσ = (1 + s) ln(1 + s) + (1 -s) ln(1 -s), s ∈ [-1, 1].
Note that F is extended by continuity on [-1, 1]. Let us introduce the set

C = { z ∈ H 1 (-L, L) : z(x) ∈ [-1, 1] for all x ∈ [-L, L] } ,
which is a closed convex subset of H 1 (-L, L). We have used the continuous Sobolev imbedding (see, e.g., [START_REF] Brezis | Analyse fonctionnelle[END_REF])

H 1 (-L, L) ⊂ C 0 ([-L, L]). ( 4.5)
It is easily seen that E K is continuous, coercive and strictly convex on C. Moreover,

E K is symmetric in the sense that E K (x → -z(-x)) = E K (z)
. Thus, we have (see, e.g. [START_REF] Ekeland | Convex analysis and variational problems[END_REF]):

Proposition 4.1. The functional E K has a unique minimizer y K on C, which is odd. If y ∈ C 2 ([-L, L]
) is a classical solution of (4.1), then y = y K .

Proof. We prove that E K is coercive, i.e. E K (z) → +∞ as ∥z∥ H 1 → +∞ with z ∈ C (the other claims are immediate). We have

f ′′ (s) = 2 1 -s 2 ≥ 2, s ∈ (-1, 1). (4.6)
Thus, f (s) ≥ 2s for s ∈ [0, 1) and

F (s) ≥ s 2 for s ∈ [-1, 1]. This shows that E k (z) ≥ 1 2 ∥z∥ 2 H 1 (-L,L) -Kz(+L) + Kz(-L).
The linear terms are controlled by the quadratic term thanks to the Sobolev imbedding (4.5). This proves that E K is coercive.

On studying the phase portrait of (4.1), the following result can be proved.

Proposition 4.2.

There exists K + > 0 such that for 0 ≤ K < K + , problem (4.1) has a unique solution y ∈ C 2 ([-L, L]) and for K > K + , problem (4.1) has no classical solution. For K = K + problem (4.1) has a unique solution

y + ∈ C 2 ((-L, L)) ∩ C 1 ([-L, L]) which satisfies y + (-L) = -1 and y + (+L) = +1.
Proof. We have already seen that any classical solution y is equal to y K . Conversely, we note that any odd solution satisfies y(0) = 0, so we consider the maximal solution of -y ′′ + f (y) = 0 (4.7)

with initial condition y(0) = 0 and y ′ (0) > 0 (the case y ′ (0) < 0 is obtained by symmetry and for y ′ (0) = 0, the solution is y ≡ 0 on (-∞, +∞)). By symmetry, any such solution is odd and defined on a maximal interval (-T + , T + ). On multiplying (4.7) by y ′ and integrating, we find that (y ′ ) 2 /2 -F (y) = C for some positive constant C. Since y ′ (0) > 0, y is increasing on [0, T + ) and y ′ = √ 2C + F (y) on [0, T + ). Thus, T + is defined by

∫ 1 0 dz √ 2C + F (z) = T + .
This shows that T + is a continuous and (strictly) decreasing function of C, with lim

C→0 + T + = +∞ and lim C→+∞ T + = 0.
Therefore, there exists a unique C L ∈ (0, +∞) such that T + = L. This defines the solution

y + ∈ C 2 ((-L, L)) ∩ C 1 ([-L, L]) with K + = y ′ + (L)
and it shows that for 0 < C < C L , the corresponding solution y is defined on [0, T + ) which contains [0, L]: we have a classical solution of (4.1). For C > C L , T + < L and we have no classical solution of (4.1).

- whose existence is asserted by Proposition 4.1. We have used a finite difference approximation with a uniform subdivision x i = -L + ih (i = 0, 1, . . . , 2M ) of [-L, L] with step size h = L/M . A discrete version of the energy E K (z) with the values (z i ) 0≤i≤2M approximating z(x i ) was used.

We have chosen the parameters L = 2, M = 20, and the minimization was performed with the Matlab ® nonlinear programming solver fmincon with constraints -1 ≤ z i ≤ 1. This solver finds the minimum of a function with inequality constraints and does not require the gradient of the function (the gradient would be singular for z i = ±1).

We have found that the critical value K + in this case is approximately K + = 1.7. For K < K + , the solution y K is regular (Proposition 4.1) and for K > K + , we have y K = y K + (Theorem 4.3).

Nonlinear solvers for the minimization of a function are known to be more efficient if a gradient is available (see, e.g., [START_REF] Chen | Positivity-preserving, energy stable numerical schemes for the Cahn-Hilliard equation with logarithmic potential[END_REF]). In view of the 2d computations, we have therefore computed numerically the solutions of the two regularized problems In this case, we replace f by the regularization

y ε K = arg min z∈H 1 (-L,L) E ε K (z) (4.8)
f ε (s) =     f (s) if |s| ≤ 1 -ε, f (1 -ε) if s > 1 -ε, f (-1 + ε) if s < -1 + ε, (4.10)
which is continuous and odd on R. We set

F ε (s) = ∫ s 0 f ε (σ)dσ, s ∈ R, (4.11)
which is of class C 1 on R and even. The energy function E ε K is defined as E K (cf. (4.4)) with F replaced by F ε . As a consequence, E ε K is of class C 1 on H 1 (-L, L). Moreover, E ε K is coercive and strictly convex on H 1 (-1, 1), so that the minimizers y ε K and ỹε K are well-defined and unique. The function y ε K solves a nonconstrained convex optimization problem, whereas ỹε K solves a constrained convex optimization problem.

For the numerical simulations presented in Fig. 2, we have set K = 3 which corresponds to a singular case for y K (cf. Fig. 1); this singular solution y K , which corresponds formally to the case ε = 0, is represented in black in Fig. 2 (left and right).

The solution y ε K is represented for several values of ε on the left of Fig. 2. It was computed by a finite difference approximation of problem (4.8) as previously, with a uniform subdivision of step size h = 2L/2M for M = 20. The Matlab ® minimization solver fminunc was used. As ε gets closer to 0, we observe that y ε K converges to y K . However, the speed of convergence is slow.

In contrast, in Fig. 2 (right), we have computed the solution ỹε K for a reasonable value ε = 0.1, and there is no visual difference between y ε K and y K . In comparison with the calculation of y K , the fmincon solver could be used efficiently for ỹε K with a much larger number of unknowns, namely 201 (corresponding to M = 100), thanks to the gradient available for the regularized discrete problem. 4.3. Interpretation of the singularity. Let K > K + . Formally, we have y ε K → y K as ε → 0. We also have y K = y K + (Theorem 4.3). This means that (y ε K ) ′ goes from K + to K on a small interval near x = +L, and at the limit ε = 0, y ′ K has a discontinuity at x = +L, which jumps from K + to K. Thus, we expect that y ′′ K has a Dirac measure at x = +L (and also at x = -L, by symmetry). We make this formal argument rigorous in Corollary 4.4.

For this purpose, we use the notion of variational solution which was proposed by Miranville and Zelik [START_REF] Miranville | The Cahn-Hilliard equation with singular potentials and dynamic boundary conditions[END_REF] in the context of the Cahn-Hilliard equation. First assume that y is a classical solution of (4.1). On replacing φ by y -φ in (4.3), we have

∫ L -L y ′ (y ′ -φ ′ )dx ∫ L -L f (y)(y -φ)dx = K[y(+L) -φ(+L)] -K[y(-L) -φ(-L)],
for all φ ∈ H 1 (-L, L). The monotonicity of f yields

∫ L -L f (y)(y -φ)dx ≥ ∫ L -L f (φ)(y -φ)dx.
We also note that

∫ L -L y ′ (y ′ -φ ′ )dx ≥ ∫ L -L φ ′ (y ′ -φ ′ )dx.
Thus, if y is a regular solution of (4.1), we have

∫ L -L φ ′ (y ′ -φ ′ )dx + ∫ L -L f (φ)(y -φ)dx ≤ K[y(+L) -φ(+L)]
-K[y(-L) -φ(-L)], (4.12)

for every test function φ.

We say that y is a variational solution to problem (4.1) if y ∈ H 1 (-L, L),

-1 < y(x) < 1 for all x ∈ (-L, L), f (y) ∈ L 1 (-L, L)
and (4.12) holds for all φ ∈ H 1 (-L, L) such that f (φ) ∈ L 1 (-L, L).

From this definition, we may deduce: Theorem 4.3. Let K + and y + be as in Proposition 4.2. Then for K ≥ K + , we have y K = y + .

Proof. Since y + is a regular solution of problem (4.1) for K = K + , it is also a variational solution, henceforth it satisfies (4.12), that is

∫ L -L φ ′ (y ′ + -φ ′ )dx + ∫ L -L f (φ)(y + -φ)dx ≤ K + [y + (+L) -φ(+L)] -K + [y + (-L) -φ(-L)], (4.13) for all φ ∈ H 1 (-L, L) such that f (φ) ∈ L 1 (-L, L). If φ is a test function, then φ is continuous on [-L, L] by (4.5) and φ(x) ∈ [-1, 1] for all x ∈ [-L, L]. Thus, we have [y + (+L) -φ(+L)] = 1 -φ(+L) ≥ 0 and -[y + (-L) -φ(-L)] = -[-1 -φ(-L)] ≥ 0.
This shows that (4.13) also holds if we replace K + by any value K ≥ K + . In other words, y + is a variational solution of problem (4.1) for all K ≥ K + . The variational solution is unique [START_REF] Cherfils | Long time behavior of the Caginalp system with singular potentials and dynamic boundary conditions[END_REF], and it can also be shown by a regularization process that the energy solution y K is also the variational solution (see Theorem 2.5). Thus, y K = y + for all K ≥ K + .

As a consequence, we have (compare with (4.3)):

Corollary 4.4. Let K ≥ K + . Then y K satisfies ∫ L -L y ′ K (x)φ ′ (x)dx -Kφ(+L) + Kφ(-L) + ⟨µ K , φ⟩ C 0 ([-L,L]) = 0, for all φ ∈ H 1 (-L, L), where µ K ∈ [C 0 ([-L, L]) ′ is the measure defined by µ K = f (y K ) + (K -K + )δ +L -(K -K + )δ -L .
Here, δ ±L denotes the Dirac measure at x = ±L. We recall that f (y K ) ∈ L 1 (-L, L).

Proof. Since y K = y + and y + is a classical solution of problem (4.1), equation ( 4

.3) shows that ∫ L -L y ′ K (x)φ ′ (x)dx -K + φ(+L) + K + φ(-L) + ∫ L -L f (y K (x))φ(x)dx = 0,
for all φ ∈ H 1 (-L, L). The claim follows.

Remark 4.5. The apparition of a Dirac measure has been analyzed for a closely related problem in [START_REF] Pierre | Global existence via a multivalued operator for an Allen-Cahn-Gurtin equation[END_REF]. It can be related to the fact that the maximal monotone operator involved in this problem is multivalued (cf.

2d numerical results

We have performed numerical simulations in two space dimension for the Caginalp system (1.1) on the rectangle Ω = [0, L x ] × [0, L y ] with L x = 8 and L y = 4, with periodic boundary conditions on the left/right sides and dynamic boundary conditions on the upper/lower sides (this type of domain is also known as a "slab").

For the space discretization, we used a finite element approximation with piecewise continuous (P 1 ) finite elements both for the order parameter u and for the temperature w. The rectangle was divided into 6400 triangles obtained by dividing the rectangle into 3200 squares of side h = L x /80 = L y /40 = 1/10, each square being itself divided into two right-angled isosceles triangles along the south-west/north-east diagonal.

In problem (1.1), we chose the logarithmic function f defined by (1.2) and the parameters were set equal to λ = 3 and α = 0. The function g was set equal to g(s) = s -β where the constant β will be specified below.

For the time discretization, we used the doubly splitting scheme (3.2)-(3.4) with a continuous regularization of f , namely the approximation f ε given by (4.10). The convex minimization problem was performed with the constraint that the solution u takes values in [-1, 1], as in (4.9). The s were made with the FreeFem++ software [START_REF] Hecht | New development in FreeFem++[END_REF]. For the optimization problem, the Truncated Newton algorithm from the ff-NLopt package was used. It requires the objective function and its gradient. The algorithm was stopped when the variation of the objective function was smaller than 1e -6. The maximal number of evaluations of the objective function was set to 200 (but it was actually never reached).

The initial values were w 0 (x, y) = 0 and u 0 (x, y) = ( 0.1 cos

( 4πx L x
) + 0.05 sin

( 6πx L x )) ( 1 - ( y L y ) 2 )
in Ω.

5.1. Numerical validation with a regular solution. For the numerical validation of our scheme, we used the parameter β = 0.8 in the function g, so that g(s) = s -0.8. With this choice, we have g(-1) = -1.8 < 0 < 0.2 = g(1) and the sign condition (cf. (2.21)) is satisfied. This guarantees that the solution is classical. We compared the solution obtained by our scheme and the solution of (1.1) computed with the linearly implicit scheme used in [START_REF] Cherfils | Long time behavior of the Caginalp system with singular potentials and dynamic boundary conditions[END_REF] (this scheme is implicit for the linear terms and explicit for the nonlinear terms). 1 shows the L 2 -error between the solution u DS from our doubly splitting scheme and the solution u LI from the linearly implicit scheme at the final time T = 0.5. The time step was chosen as dt = T /(10 × 2 m ) with m ∈ {0, 1, . . . , 5}. Since both schemes are first order in time, the difference ∥u DS (T ) -u LI (T )∥ L 2 (Ω) is expected to be first order at most. The ratio of consecutive errors between two consecutive time steps, which is very close to 2, is consistent with the first order approximation. We point out that for a regular solution, the linearly implicit scheme is much faster than our doubly splitting scheme. This is illustrated in Table 2, where the CPU computational time is represented. We used a laptop for which the computation of u LI (T ) with the time step dt = T /10 takes about 6 seconds. Since the CPU time depends on the computer, in Table 2, the CPU time is normalized by comparing every computation to this case (for which the CPU time is therefore set to 1). For the linearly implicit scheme, when the time step is divided by 2, the CPU time is multiplied by 2, as expected. For the doubly splitting scheme, the most greedy part is the minimization algorithm. The computational time could be reduced by using a second order method involving the hessian of the objective function and not only the gradient (in such a case, a C 1 regularization of f would be used, yielding a C 2 regularization of F ). step is dt = 0.01. We have g(-1) = -4 < g(1) = -2 < 0 so the sign condition (2.21) is not satisfied. The solution u reaches the value +1 somewhere on the boundary at time t = 0.71 and is singular for t ≥ 0.71. The isovalues of u are represented in Figures 3-6 at times t = 0, t = 0.10, t = 0.71 and t = 5.00 (in these Figures, the Matlab ® software was used for the visualization). We observe that the solution converges to a singular steady state, which is constant along the x-direction.

In order to visualize the evolution in the y-direction, we have represented in Figure 7 the solution u on the segment {x = 2} × [0, L y ] as time grows. For this particular section of u, the value +1 is reached at time t = 0.73. We observe that, afterwards, u continues to evolve in the bulk, until it reaches a singular steady state.

Figure 1 . 4 . 2 .

 142 Figure 1. Solution y K for different values of K 4.2. Numerical resolution. In Figure 1, we have computed numerically the solution y K = arg min z∈C

Figure 2 .

 2 Figure 2. Solutions of the regularized problem without constraint (y ε K , left) and with constraint (ỹ ε K , right)

  (2.8)).Corollary 4.4 shows that for problem (4.1), the support of the singular part belongs to the boundary of the domain and corresponds to a discontinuity of the normal derivative of the solution. The latter phenomenon has been described for the evolutionary 3d Caginalp problem in terms of the variational solution of the problem in[START_REF] Cherfils | Long time behavior of the Caginalp system with singular potentials and dynamic boundary conditions[END_REF] Theorem 5.2].

Figure 3 .

 3 Figure 3. Initial condition u 0

Figure 4 .

 4 Figure 4. Solution u(t) at time t = 0.10

Figure 5 .

 5 Figure 5. Solution u(t) at singular time t = 0.71

Figure 6 .Figure 7 .

 67 Figure 6. Stationary solution (u(t) at time t = 5.00)

Table 1 .

 1 L 2 -error and ratio of consecutive errors vs time step Table

	m (cf. time step)	0	1	2	3	4	5
	L 2 -error	0.0240 0.0124 0.0063 0.0032 0.0016 0.0008
	ratio	1.94	1.97	1.97	2	2

Table 2 .

 2 Normalized CPU time vs time step for the linearly implicit (LI) scheme and the doubly splitting (DS) scheme

	m (cf. time step) 0	1	2	3	4	5
	LI scheme	1	2	4	8	16 32
	DS scheme	165 262 305 381 511 489

Acknowledgements

This paper was initiated while the first author was visiting the "Laboratoire de Mathématiques et Applications" in Poitiers, whose hospitality is gratefully acknowledged. The visit of the first author was supported by the AUF (Agence Universitaire de la Francophonie) in the framework of the project "Etude des modèles mathématiques de pollution des eaux de surface et souterraines".


François-Xavier Coudert
email: fx.coudert@chimie-paristech.fr

Reproducible research in computational chemistry of materials

The reproducibility of experimental findings is a crucial underlying tenet of the scientific method, and the awareness of its central place in the advancement of knowledge has been long ) in the 11 th century and later to the scientists of the Renaissance. While the way we perform and report research has evolved immensely, reproducibility remains as vital today as it was then -and maybe more, in a world where 2.5 million scientific articles are published worldwide every year. A quantitative analysis of the causes of retraction of research articles, published in 2012, [START_REF] Grieneisen | A Comprehensive Survey of Retracted Articles from the Scholarly Literature[END_REF] concludes that 22% of retractions cite a "known artifact" in the data or "unexplained irreproducibility" as the cause of the retraction, without any research misconduct or ethical violation. This statistic highlights the need to improve the current standards for the reproducibility of published works, by creating and promoting policies for reproducible science. This parallels the recent advances by journals, professional associations, academic institutions, and funding agencies in the domain of publication ethics (see, for example, the guidelines developed the Committee on Publication Ethics 2 and nowadays adopted by many scientific journals).

In this landscape the swift rise of computational science, which has led to exciting developments and breakthroughs in all fields of academic research, requires us to reinterpret the traditional understanding of reproducibility to apply to computational results from in silico experiments.

Computer calculations being a finite series of operations on fully deterministic machine, they are in theory fully replicable, and can be rerun to yield bit-for-bit identical results. However, in practice this does not translate at all into full reproducibility of computations, despite the modest cost of online hosting data, owing to the complexity of the modern hardware and software stacks. I have, like many colleagues, been frustrated a few times in the past by difficulties in reproducing results from the published literature -or by obtaining structures or data related to published findings, for comparison with my own work. This is the case even when the authors would be fully willing to share their data at the time of publication: over time the data might have been lost, or the person working on it gone, and there is no real incentive to go dig up files from 5 or 10 years backespecially in cases where the group does not work in this area anymore.

Browsing through a recent issue of a flagship materials chemistry journal gives a good idea of how important computational methods have become in our field, and what progress we still have to make towards reporting these studies in a way that promotes reproducibility. Out of 41 full research articles in an issue, I counted 7 that were exclusively or predominantly theoretical in nature (17%), with 13 more that included some computational results in figures or tables (thus a total of 49%). For most of the papers, there was little to provide any help to a researcher willing to reproduce the calculations: the molecular and crystal structures discussed were not provided except as snapshots in figures; input files for calculations were not provided. In a few cases, where computational work was one of several characterization methods, the level of description of the work performed (whether in the main text or as supporting information) actually appeared too low to allow reproducing the calculations. Three papers provided full structural models, three others reported extensive datasets of computed properties, and one gave the source code used for analysis.

However, in both cases, these were available only in PDF format over multiple pages.

Here, I summarize the different aspects of research reproducibility in computational chemistry and materials science, and delineate the current practices of the community (both researchers and publishers) in this respect.

Opening the data

The first necessary step toward research reproducibility and the collective advancement of knowledge is that scientists make available the results of their experiments and calculations.

Numerical quantities measured are presented in the text, in tables, or in graphs. Moreover, all structures -molecular or crystallographic -discussed should be made available, not only in graphical form (e.g. as 3D representation) but as structure files. Most journals nowadays include this requirement in their instructions to authors, although it is unfortunately not always enforced.

It is also good practice for the structures to be given as machine-readable files in standard formats (XYZ, PDB, or CIF files) rather than to list coordinates in a PDF document as supporting information. This lowers the barrier for others to visualize the data or reuse it in future works, as well as improves the discoverability and indexing by future data-mining projects. Some computational studies can produce a large amount of data (chemical structures as well as computed properties), e.g. in works of high-throughput screening or large scale calculations on existing structure databases. This requires thorough curation of the data to produce a usable database, including all properties produced in a standard, machine readable format (such as XML or JSON) accompanied with the documentation of the format itself. Care should also be taken to make sure that the database published contains all relevant metadata, for example the calculation methods and method parameters employed (if they vary between structures) or the conditions of the calculation (elapsed time, platform on which it was run, etc.). When possible, this data can then be cross-linked with the earlier database(s) from which it was build. An excellent example of this practice is the recent inclusion of computed elastic properties for inorganic crystalline compounds [START_REF] De Jong | Charting the complete elastic properties of inorganic crystalline compounds[END_REF] from the Materials Project: [START_REF] Jain | The Materials Project: A materials genome approach to accelerating materials innovation[END_REF] the computed dataset was first published online (DOI: 10.5061/dryad.h505v) and linked to in the paper, and then integrated into the Materials Project web-based explorer. It original contained computed data for 1,181 materials, while that number has since then gone up to 4,375.

Finally, let's note that the choices made for hosting the data online impact its long-term availability. While hosting data on one's group website may seem practical, it does not guarantee availability in the long term. The use of larger-scale institutional repositories should be preferred -or not-for-profit repositories with reasonable guarantees of long-term storage and independent archiving. Getting a stable URL, that will not change over time, is also important: some online services, such as Zenodo, 5 offer Digital Object Identifiers (DOI) for data stored in public GitHub repositories. [START_REF]7 Dealing with data[END_REF] Recent years have seen a large improvement in data availability, coming both from individual author practices but also from changes in several journals' policies, which now request statements declaring the accessibility of the data and its location. 7 Statements that "data is available upon request" may rapidly become a thing of the past -only to be used in special circumstances the data cannot be freely released.

Opening the input and output files

In addition to the publication of curated data that is presented and discussed in the paper, full reproducibility of computational studies requires to make available to the research community at large the input files that were used to start the computations. A complete set of input filesincluding all parameters and initial configuration of the system -contains all information about a given computation, includes many "technical" parameters that may not always be included in an article's methodology section. This is useful to both the reader and the reviewer, who is tasked with judging whether the methods applied are sound. Moreover, it provides information directly in machine-readable form, lowering the bar for replication studies and avoiding potential human errors -avoiding for example to rekey an entire force field input file from article tables.

In addition to input files, in a context where studies routinely use large amounts of highperformance computing (HPC) resources and storage is relatively cheap, the next logical step in improving reproducibility is to store and publish the output files (or raw data) from the computations -in all or part. This allows readers to study the details of the computations ran, even if they do not have access or do not wish to spend large CPU time rerunning the calculation.

It also enables other researchers to use data from the raw output that may not have been exploited by the original authors, but could be of interest for other purposes -atomic charges for analysis or force field derivation, vibrational eigenvectors, etc. This usefulness is, however, balanced by the fact that raw calculation outputs can take up large amounts of disk space: neither publishers nor authors would be willing pay the cost of long-term hosting for unreasonably large files, and the output files thus need to be curated (trimmed or selected) before being published. We also note in passing that some restrictive software licenses do not allow you to publish output filesrequiring you, for example, to redact execution times or computational performance data.

A novel project has been born, whose goal is to build a distributed database of computational chemistry results. ioChem-BD 8 has a two-pronged focus on converting data from several codes into a common data standard, as well as allowing data management, search and manipulationall the way to publishable supporting information files. It allows storage of the data, a graphic user interface to manipulate it, and APIs to connect it to other databases.

Opening the software

Finally, we turn our focus now to the core of the computational chemistry research: the software.

Indeed, even with the publication of computation input files, full reproducibility requires the use of the same software. Even the same method, implemented in two different codes, can lead to different results. This means that the use of commercial software or unpublished in-house codes can restrict the reproducibility of a computational study by the research community. I consider this an encouragement to promote and use software that is either open source or at least freely available to the academic community. Furthermore, in the case were one uses commercial (or not publicly available) software, there is an even stronger incentive to publish online the raw outputs of the computations, thus mitigating the issue.

The remarks above are true not only of large computational software packages, such as quantum chemistry packages or molecular simulation suites, but also of the smaller pieces of software created during the course of the research: analysis codes, post-processing tools, visualization scripts, etc. Those are sometimes overlooked, yet often crucial results depend on their correctness -and the details of analysis methods are not always fully documented in the articles themselves.

Finally, let us note that even for software packages whose source is freely available, differences in the exact results obtained can vary depending on software version, but also on the nature of the hardware used, on the other software used on the computer (operating system, compiler, mathematical library, etc.), as well as on runtime parameters such as number of CPU cores or HPC nodes used. Thus reproducibility of research and publication of data is particularly important in computational sciences, [START_REF] Lejaeghere | Reproducibility in density functional theory calculations of solids[END_REF] even when exact replicability of the simulations is not technically possible -because the calculations used hardware that has become obsolete, for example.

In conclusion, the ease and low cost with which we can nowadays host reasonable amounts of data online is an encouragement to publish more data accompanying scientific articles, going beyond the traditional "supporting information" files to include input files, output files of computations as well as source code. This will improve reproducibility of existing work and strengthen confidence of the wider public in the role of computational chemistry tools in materials science.

data input files raw output source code


François-Xavier Coudert

Correcting the scientific record: retraction practices in chemistry and material sciences

currently form the cornerstone of the modern scholarly publication system, and guarantee the dissemination of research findings through the worldwide, ever-increasing community of researchers. Collectively these published works, stamped with the seal of approval of a review by the authors' peers, form the scientific record -

Despite the healthy role of retractions in preserving the scientific record, and while erroneous data can be the result of a good faith mistake, there is definitely a stigma associated with the retraction of a paper. COPE guidelines state that "The main purpose of retractions is to correct the literature and ensure its integrity rather than to punish authors who misbehave", 3 but previous work has shown a notable resistance to admitting error in scientific papers. 4 The term retraction is too often associated with research misconduct, giving it a negative connotation for authors. [START_REF] Barbour | Amending published articles: time to rethink retractions and corrections? F1000Research[END_REF] This is particularly true in a highly competitive environment, where academics are driven to publish often and produce high-impact papers: Jin et al. showed that retractions have a negative effect on citations for earlycareer researchers. 6,[START_REF] Lu | The Retraction Penalty: Evidence from the Web of Science[END_REF] The same argument can also be made for the publishers, who may fear a dent in the reputation of the journal. Thus, none of the actors involved have any direct incentive to retract a paper.

In this context, and while the number of retractions is rising, [START_REF] Cokol | Retraction rates are on the rise[END_REF][START_REF] Steen | Why Has the Number of Scientific Retractions Increased?[END_REF] there is relatively little information available about retractions and retracted papers, beyond the retraction notices infrequently published by journals. There is no central repository or authoritative database that can be easily queried -although the Retraction Watch website, which covers the topic of retractions and publication ethics in general, has been collating such a database. 10 Previous systematic studies have focused on retractions in specific fields, and in particular in medicine [START_REF] Wager | Why and how do journals retract articles? An analysis of Medline retractions 1988-2008[END_REF][START_REF] Redman | Empirical developments in retraction[END_REF][START_REF] Moylan | Why articles are retracted: a retrospective cross-sectional study of retraction notices at BioMed Central[END_REF] -with the notable exception of a study by Grieneisen et al. that spanned several fields of research. [START_REF] Grieneisen | A Comprehensive Survey of Retracted Articles from the Scholarly Literature[END_REF] In order to better understand the existing practices for article retractions in the chemical sciences, I have performed a systematic survey of 331 papers retracted in 2017 and 2018 and their retraction notices, publicly available on the journals' websites. This article looks at the statistics of retractions, their distribution per country, the occurrence of multiple retractions. I also provide a classification of the reasons behind the retractions, and the distribution of their occurrence.

How many papers are retracted?

A search of the Elsevier's Scopus database 15 for retraction notices, in the fields of chemistry, materials sciences, and chemical engineering, resulted in a list of 331 papers that were retracted during the years 2017 and 2018 (see the Methods section for details). This is a very small number compared to a total of 1,114,476 papers published in those fields in the same two years, with 3 papers retracted for 10,000 papers published.

The geographic distribution of the retracted papers is given in Figure 1. The distribution generally reflects the number of papers published by each country. We see that 31.4% of retracted papers involve a Chinese affiliation, which is in line with the fact that 30.9% of papers published in 2017-2018 have a Chinese co-author. The same is true of the United States, with 14.2% of retractions (15.1% of published papers), South Korea (4.5% vs. 4.0%), and other countries lower in the listalthough with smaller absolute numbers, the statistics become noisier. The ratio of retraction in China is of particular interest here, given the country's cash-per-publication reward policy 16 has been suspected to drive researchers toward, at best, the production of more papers of lower scientific significance, or at worst a negative impact on academic ethics. [START_REF] Fuyuno | Cash for papers: putting a premium on publication[END_REF] The data show that this does not, for the field and period studied, lead to an increase in the ratio of retraction.

There are, however, countries whose "share" of retractions exceed their share of published papers. Indian affiliations account for 17.2% of retractions in our database, compared to 7.3% of published papers worldwide -a difference of a factor of 2.4. But the country featuring for the largest rate of retraction is Iran, which is involved in 11.2% of retractions compared to 2.7% of published papers. The ratio of retracted papers from Iran is 4.1 times larger than from the rest of the world. I note that while a small number of authors contribute heavily to those number, they do not -in themselves -appear to account entirely for the higher rate. Analysis over a larger time period would be necessary to confirm these trends.

Figure 1. Distribution of author affiliations, per country, for retracted papers. The histogram lists

the percentage of retracted papers where at least one author is from a given country. Total of percentages exceeds 100% because papers can be written by co-authors from different countries.

How long after their original publication are papers retracted? The distribution of retraction time, presented in Figure 2, is well-spread. 10% of papers are retracted within 3 months of the original publication. In several cases the issues were discovered even before papers were "formally" published as part of a journal issue, while the accepted versions of the manuscripts were posted online. On the other hand, many papers take longer to retract: the median time to retraction in this dataset is 24 months. But the distribution has a long tail, the tardiest retraction observed is for a paper published in June 2000, and retracted 18 years later, in February 2018 (DOI: 10.1007/s12540-018-0048-0).

Figure 2. Distribution of time to retraction, in months.

Looking at the distribution in citations of the retracted papers (Figure 3), we see also a wide variety. This is linked, to a large extent, with the large distribution of time to retraction. Since we look at papers retracted in 2017 and 2018, the vast majority of citations reported were made before the retraction -it would be interesting, in future work, to look at older retractions to see if those papers are still cited by the community. The median number of citations is quite low: 4, while the average is 10.7. The most cited paper in the dataset was cited 121 times in the course of 13 years (DOI:

How common are multiple retractions?

The past few years have seen the revelation of a certain number of cases of high-profile scientific misconduct, where the discovery of unethical practices by an author or research group leads to retraction of several papers. From the 2017-2018 retractions analyzed here, we clearly identify some cases of large-scale retractions, all of which have been reported in specialized media. The largest "clusters" of retractions in that two-year period include 18 retractions for R. Madhuri and P. K. Sharma (Indian Institute of Technology, Indian School of Mines in Dhanbad), [START_REF] Oransky | Pair of nanotech researchers up to at least two dozen retractions[END_REF] 10 for G. Tang and X. Ye (Tsinghua University in Beijing), [START_REF] Marcus | Graduate student in China stripped of PhD after investigation that led to a dozen retractions[END_REF] 9 for A. S. Elahi and coworkers (Islamic Azad University in Tehran), [START_REF] Stern | After Elsevier knew an author faked reviews, it kept accepting his papers for more than a year[END_REF] 7 for S. Thakur and R. B. Tokas (Bhabha Atomic Research Centre in Mumbai), [START_REF] Stern | [END_REF] and 7 also for S.-G. Shin (Kangwon National University in Chuncheon). [START_REF] Stern | [END_REF] These retractions were, in every case, linked to ethical shortcomings.

These cases of large-scale retractions are well publicized, and they account for a nonnegligible fraction of the retracted papers. We count 88 retracted papers involved in multiple retractions (3 or more papers retracted with at least one common co-author), which makes for 27% of the retractions. On the other hand, out of 1093 authors of retracted papers in 2017 and 2018, 92% (1003 authors) are involved in a single retracted paper, and only 31 authors (3%) are involved in the 27% of multiple retractions.

Why are papers retracted?

To see why papers are being retracted, I have classified the reasons being listed in each retraction notice (see Methods for details) into four main categories. The results are shown in Figure 4, and the data is available as supporting information -including, for each paper, the full retraction notice and my classification into categories. The most frequent reason behind paper retractions is plagiarism, involved in 42% of cases. It can take the form of duplicate publication (15% of cases; the direct reuse of an entire article, sometimes after translation), self-plagiarism (45%; the reuse of significant parts of text without proper citation), or plagiarism of previously published works by other authors (40%). We note that, although nowadays publishers can now check for plagiarism with automated software tools upon submission of manuscripts, [START_REF] Butler | Journals step up plagiarism policing[END_REF] that was not always the case, and some papers retracted in 2017 and 2018 were published before this software was available. We cannot determine, from the published retraction notices, what is the main motivation of authors in plagiarizing, but previous studies have pointed out to possible contributing factors: lack of academic infrastructure and mentoring; [START_REF] Errami | A tale of two citations[END_REF] incentives based on the number of publications, whether in the form of cash bonuses [START_REF] Fuyuno | Cash for papers: putting a premium on publication[END_REF] or career advancement. [START_REF] Citron | Patterns of text reuse in a scientific corpus[END_REF] The second largest root cause of retraction is the presence of ethical issues with the data reported, with 90 occurrences (27%). The terminology used by retraction notices -which can be written either by the publisher, the editor, or the authors -to discuss this data can be very varied. Some notices are very specific, providing a clear characterization of the issues: "research misconduct", "data fabrication", "data falsification", "inappropriately and extensively digitally manipulated", etc. Other notices are much less clear-cut, mentioning unspecific "problems", "significant problems", "extraordinary similarities", "incorrect data", "doubts cast on the data", … I have therefore broken down the "data issues" category into three different subclasses. In a majority of 60% of the cases, the data is falsified, i.e., the unethical handling of data is clearly presented. This includes cases of fabrication, invalid duplication, digital editing, etc. and accounts overall for 16% of retracted papers. In a smaller number of cases (31% of data issues, 8% of retractions), the issues with the data are not clearly explained, and the unethical nature of the publication is merely hinted at through the language of the notice. Although this is not the majority of cases, it still appears quite high -such cases should be totally avoided, to provide maximum information to the scientific community, and to obey COPE guidelines: 3 "retraction notices should always mention the reason(s) for retraction to distinguish honest error from misconduct". Finally, there were a few (8) cases where the cause of retraction is due to lack of data ownership, lack of authorization to publish the data, or disclosure of confidential or proprietary information.

Lack of authorization

Unclear authorship

Affiliation error The third cause of paper retraction is related to authorship issues: this is mentioned as a factor in 19% of retraction notices. In most cases, the issue is a lack of authorization to publish from one of the authors. In most cases, it is a co-author who has not been informed of the submission and publication of the paper -and quite often a senior co-author, supervisor or advisor, according to the retraction notices surveyed. In a few cases, this can mean the inclusion on the author list of a co-author who did not participate in the work at all, possibly adding a well-known name to the author list to improve a manuscript's chances of success. Authorship disputes leading to retraction can also take the form of a missing co-author (13 cases), i.e. someone whose contribution to the work met the criteria for authorship, yet was not added to the paper's author list -or cases where the contribution of different individuals is unclear or disputed (12 cases). In particular, there are signs in some retraction notices of authorship being sold, after review, a known practice in some countries: [START_REF] Hvistendahl | China's Publication Bazaar[END_REF] "the paper was submitted to the journal by the first author, who was not able to provide a reasonable explanation for the addition of the other seven authors to the revised version of the paper". Finally, some papers were retracted due to the presence of an incorrect affiliation (7 cases), although this is often a contributing factor but not the sole reason for the retraction. 27 It may be surprising to note that the least-common category of retractions is that of papers that were identified, after publication, to contain "honest errors", i.e. errors that are not attributed to misconduct. This represents 16% of retractions in this corpus, and they are almost evenly split between: (i) theoretical errors, in mathematical derivations or during data analysis; (ii) experimental errors, including experiments performed in invalid conditions or apparatus used in an unreliable manner; (iii) results that could not be reproduced, by the authors or other researchers, in later works. It is interesting to note that these retractions tend to be published at the request of the authors, and in many cases feature long retraction notices that fully detail the nature of the errors made (see, e.g., DOI: 10.1021/acs.macromol.7b00492 and DOI: 10.1038/nchem.2885), and the investigation. In some cases, these investigations into the reproducibility of the work have even been published, separately of the retraction, as a scientific article in their own right (see DOI: 10.1103/PhysRevB.98.079902 and Ref. 28). This really represents a "positive" type of retraction, that is to be encouraged for it benefits greatly the community.

Missing author

Self

Other features of retraction notices

In addition to these four main categories of retraction reasons, we note here a few other lessfrequently mentioned retraction reasons. The first is the occurrence of a technical error at the publisher (5 cases, 1.5%). For example, where a paper was accepted without review by mistake, or was published twice due to an error in the production chain. The second is the retraction of papers whose peer-review was insincere: I list this separately here, because while it affects 24 papers (7.2%), it is actually not common, and happened at a single publisher (Elsevier). An inquiry revealed that these papers had been accepted after a biased or fake peer-review: authors had submitted potential reviewers' contact details (including email), yet those emails accounts were not associated to the listed individuals, but manipulated (either by the authors or third parties). 29 Most publishers have now put in place a verification mechanism, to ensure the authenticity of reviews.

In addition to the reasons stated as causes for retracting papers, the notices can feature additional information about the process of the retraction. For example, some notices state who requested the retraction (author or authors, editor, publisher) -but since not all notices do, it was impossible to gather good statistics on this. I have also noticed that many notices (23 out of 331) mention the publisher being unable to contact some or all of the authors. This occurs in particular in cases where the retraction is due to ethical misconduct, probably because authors either try to ignore the issue, hoping it will go away, or may have judged that silence is the best defense.

Two other significant features found in some (but not all) retraction notices are the discussion of the individual responsibility of the authors, and the mention of possible links to institutional investigations, either pending or completed. Inquiries by academic institutions, usually the authors' principal affiliation, is mentioned in 19 notices (5.7%). It is explicitly mentioned, in several cases, that the article is retracted following an institutional inquiry into research misconduct. The reverse is quite rare: it does not appear that journals or publishers routinely inform the authors' academic institutions of the retraction proceedings.

Finally, 13% of the retraction notices studied discuss explicitly the individual responsibility of the authors of a retracted paper. This is rarely a systematic discussion of all co-authors, but more often done to exonerate one specific author ("X has stated that she was unaware of the submission") or to assign the blame to a single author ("sole responsibility lies with the corresponding author", "the corresponding author takes full responsibility and apologizes for the academic misconduct"). In other cases, such statements are used to indicate which authors endorse the retraction or not ("X and Y oppose the retraction"). This analysis of the retraction factors is in generally in line with the findings reported by Grieneisen et al. across several fields [START_REF] Grieneisen | A Comprehensive Survey of Retracted Articles from the Scholarly Literature[END_REF] (i.e., on all journals indexed by Web of Science).

Level of detail in retraction notices

The content of retraction notices can vary widely -in both size and information provided. We found at least one journal (Journal of the American Oli Chemists' Society) that appears not to publish freely the text of retraction notices (DOI: 10.1080/10837450.2018.1451433) (DOI: s11746-011-1789-9), a practice contrary to the COPE Guidelines. Other times, retraction notices contain no actual information as to the motivating factors. As an example, one retraction states: A likely contributing factor for such short notices is the pressure on journals and their editors to avoid using clear or harsh retraction language, for fear of being proven wrong in a further inquiry -or fear of legal retribution, as some authors have been known to contest retractions in law courts.

At the other end of the scale, some notices are extremely detailed, and provide a detailed analysis of the reasons behind the retraction. Many retraction notices for plagiarism provide a clear picture of what material was copied from the authors, and what were the original sources. Unsurprisingly, the notices of papers retracted for good-faith mistakes are often the most detailed, providing good insight to readers on the problems with the published work. For example, Lu et al. (DOI: 10.1021/acs.macromol.7b00492) retracted four papers due to incorrectly-performed molecular simulations, and gave a detailed account of the shortcomings of their simulation procedure in a fullpage retraction notice. Another interesting inquiry, related to a retracted paper, is the work of Ryan et al., who repeated the chemical reactions of a retracted organic chemistry paper (DOI: 10.1021/acs.orglett.7b02161) and provided a full analysis of the shortcomings in the analysis of the original paper. [START_REF] Ryan | Repetition of chemistry from a recently retracted paper. A cautionary note[END_REF] Some retraction notices even acknowledge whistleblowers (when they did not wish to remain anonymous, we suppose!) "for the considerable effort put into collecting the evidence and reporting this case of multiple publication", which I see as an encouraging note.

Conclusions

The aim of this study was to provide an overview of the retraction practices in chemistry and materials sciences. By analyzing retraction notices for a two-year period (2017-2018), I provide quantitative data on a phenomenon that is relatively common for the field as a whole -although quite extraordinary at the individual level. The repartition across all countries is relatively homogeneous, with two countries having significantly higher retraction rates: India and Iran. The causes of retraction are mainly related to misconduct, with plagiarism accounting for 42% of retractions, and data issues for 27%. Authorship issues are also common, and in particular papers being submitted without the knowledge and authorization of all authors, often by younger co-authors: this probably indicates that better education or training is needed in the workings of the academic publishing system and its rules.

While this analysis only concerns one academic field, nothing in the trends observed appears to be particularly tied to the specificities of chemistry or materials sciences. Its findings are thus probably generalizable to neighboring fields in experimental sciences, including physics and biology -although this cannot be directly confirmed, as I could not identify any systematic study on retraction practices in other fields, except medicine. There, a previous study 11 also identified plagiarism as the largest factor leading to retraction, but did not identify conflicts of authorship and lack of authorization by co-authors as a significant factor. This is difficult to interpret, though, and may indicate either that co-authorship is more formalized in the medical field, or that editors are more likely to correct a paper following authorship disputes, rather than retract it. This quantitative analysis of the recently retracted papers is clearly only a first step in providing insights into the practices of retraction in chemistry and materials sciences. In order to better understand the current status of retractions in the publishing model, it will have to be complemented with sociological investigation into the views of all parties involved, the motivations behind these views, and surveys of academics in their various roles in the chain of publication: as authors, reviewers, editors, and publishers.

Methods

I gathered data by a Scopus search on "retraction" or "retracted", for the years 2017 and 2018, limited to the fields of chemistry, materials sciences, and chemical engineering (as identified by the Scopus "subject areas"). I curated by hand the 362 results obtained, to remove the use of "retraction" in other contexts (e.g., retraction of polymer chains), and the result was a list of 331 unique retracted papers. For each paper, I gathered using the Scopus API the DOI of the retraction notice, the DOI of the original paper (which are sometimes the same), the authors' names, the affiliations (and in particular the country of each affiliation).

Because there is no simple API for the retrieval of the text of retraction notices, I retrieved them manually from each article's web page. Once the whole dataset was assembled, each retraction notice was manually classified to identify the following possible reasons stated as being the cause(s) for, or factor(s) in, the retraction: § Authorship issues:

-Lack of approval to publish from some of the co-authors -Author missing from the author list -Error in the affiliations listed -Authorship disputes or unclear authorship § Errors in the science reported:

-Experimental error -Error in the analysis, interpretation, or mathematical derivations -Work that was not reproducible by the authors or other researchers § Issues with the data -Data falsification or mishandling, clearly identified as such -"Problems" or "issues" with the data, for unclear reasons -Issues with data ownership, or authorization to publish (confidentiality, etc.) § Plagiarism -Duplicate publication, identical or near-identical papers -Self-plagiarism (unacceptable reuse of material previously published by some of the authors) -Plagiarism of other sources (not by the authors) § Publisher error § Issues of copyright § Problems with the integrity of the review process I also took note of the occurrence, in the retraction notices, of the following features: § Authors that could not be contacted or did not respond during the retraction process § Mention of an institutional inquiry into research misconduct by one or more of the authors' affiliations § Statement detailing the individual responsibility or role of the authors in the publication of the paper, or during the retraction process.

All data gathered, including the classification of each retracted article, is available as supporting information and in the repository at https://github.com/fxcoudert/citable-data/

  /j.biortech.2017.08.105)… and was retracted as a duplicate of another publication by the same author (which itself has gathered 130 citations).

Figure 3 .

 3 Figure 3. Distribution of citations of retracted papers.

Figure 4 .

 4 Figure 4. Breakdown of the reasons listed for retraction, classified from the text of the retraction notices. The sum of the four main categories exceeds the total number of papers, because retractions can list several contributing factors in the publisher's decision.

"

  We, the Editor and Publisher of Mineral Processing and Extractive Metallurgy Review are removing the following article: [article title]" (DOI: 10.1080/10837450.2018.1451433) Another manages to give slightly more information while keeping readers in the dark: "This article has been retracted at the request of the authors because of a business decision by their employer." (DOI: 10.1016/j.jcrysgro.2018.09.045).

Table 1 .

 1 Breakdown of the reasons listed for retraction, classified from the text of the retraction notices.

	Plagiarism	139
	including: Duplicate publication	21
	Self-plagiarism	62
	Plagiarism of others	56
	Data issues	90
	including: "Problematic" data	28
	Falsified data	54
	Ownership or disclosure issues	8
	Authorship issues	62
	including: Missing authorization to publish	30
	Missing author	13
	Unclear authorship	12
	Affiliation error	7
	Honest errors	54
	including: Error on the theory or data analysis	21
	Experimental error	17
	Work not reproducible	16
	Other factors mentioned in retraction notices
	Publisher error	5
	Copyright issue	3
	Abuse during review	24
	Unable to contact one or more co-authors	23
	Institutional inquiry

Acknowledgments

Nicolas Jonckheere is gratefully acknowledged for critical reading and comments on a draft of this paper.

Conflicts of interest

FXC is Associate Editor for Adsorption Science & Technology (SAGE Publications), and Editorial Board Member for Communications Chemistry (Springer Nature). He has not, in these capacities, handled any retraction procedure.