Computer-assisted equilibrium validation for the diblock copolymer model

The diblock copolymer model is a fourth-order parabolic partial differential equation which models phase separation with fine structure. The equation is a gradient flow with respect to an extension of the standard van der Waals free energy functional which involves nonlocal interactions. Thus, the long-term dynamical behavior of the diblock copolymer model is described by its finite-dimensional attractor. However, even on one-dimensional domains, the full structure of the underlying equilibrium set is not fully understood. In this paper, we develop a rigorous computational approach for establishing the existence of equilibrium solutions of the diblock copolymer model. We consider both individual solutions, as well as pieces of solution branches in a parameter-dependent situation. The results are presented for the case of one-dimensional domains, and can easily be implemented using standard interval arithmetic packages.


1.
Introduction. Over the years, mathematical models for phase separation phenomena and the resulting questions of pattern formation have proved to be fertile ground for the development of new mathematical tools and techniques, and in many cases these developments have provided important new insight into associated questions in materials science [13,14,15]. One of the particularly interesting directions has been the discussion of models which are not completely local in nature, but include in one way or another nonlocal effects [6,16,17]. While in some models the nonlocal effects show up directly in the evolution equation, in other models they enter through an associated energy. One of the latter situations occurs in the study of diblock copolymers, which are formed by the chemical reaction of two linear polymers, or blocks, which contain different monomers. These blocks are often thermodynamically incompatible and therefore tend to separate. Yet, if these blocks are already covalently bonded, they cannot separate on a macroscopic scale without adopting an entropically unfavorable extended configuration -and this leads to the phenomenon of microphase separation, where the two block types separate on a mesoscopic scale.
One of the basic models for microphase separation in diblock copolymers has been described in [33], and in its original form was proposed by Ohta and Kawasaki [34] and Bahiana and Oono [2], see also [7,8]. If we consider a material constrained to the bounded domain Ω ⊂ R d , then their model considers a free-energy functional for the relative macroscopic monomer concentration u, i.e., the difference between the two monomer volume fractions, which is given by In this formula, the real number µ = Ω u dx/|Ω| denotes the average mass of u, and the nonlinear function W is a double-well potential with global minima at ±1.
In the following, we always consider the standard double-well potential W (u) = (u 2 − 1) 2 /4. Notice that the energy functional (1) is just the standard van der Waals free energy, but with an additional nonlocal term which involves the square root of the inverse Laplacian on a space with zero total mass. While there are different ways in which one can associate dynamics to (1), the diblock copolymer model in its standard form considers the evolution equation which is based on the gradient in the H −1 -topology. In this formulation, the value of u describes the local material composition in the following way: Values of u close to +1 are interpreted as only component A being present at a point x ∈ Ω, and the value −1 indicates that only component B is present; values in between correspond to mixtures of the two components, with zero representing an equal mixture. The parameter µ denotes the average mass of the mixture, and the two remaining parameters λ > 0 and σ ≥ 0 are dimensionless parameters, which are described in more detail in [22]. We would like to point out that for σ = 0, the diblock copolymer model (2) reduces to the well-known Cahn-Hilliard equation, which serves as a basic model for the phase separation phenomena spinodal decomposition [26,27,42,43,45] and nucleation [4,5,11], see the survey in [14]. As a dissipative parabolic equation, the long-term dynamics of the diblock copolymer model (2) is described by its associated global attractor. Furthermore, since (2) is a gradient system with respect to the energy (1), the global attractor consists of equilibrium solutions, as well as heteroclinic solutions between them. Thus, at least in principle one should be able to gain significant insight into the long-term dynamics by understanding the structure of the equilibrium set -and since the diblock copolymer model arises from the Cahn-Hilliard model as a regular perturbation, one might be tempted to assume that their attractor structure closely resembles one another, at least for small values of σ. Yet, nothing could be further from the truth.
Consider for example the case of the one-dimensional domain Ω = (0, 1). For the Cahn-Hilliard model the complete equilibrium set is known [21], and it shows for example that for mass µ = 0 and sufficiently large λ > 0 there are exactly two global energy minimizers, which globally attract almost all solutions of the equation. In contrast, in the diblock copolymer case σ > 0, the number of stable equilibrium  (2) on the domain Ω = (0, 1) and with σ = 6. The left diagram is for total mass µ = 0, while the right diagram is for µ = 0.3. The solution measure for both diagrams is the L 2 (0, 1)-norm of the solutions. The solution branches are colorcoded by the numerically determined Morse index of the solutions, and the colors black, red, blue, green, magenta, and cyan correspond to indices 0, 1, 2, 3, 4, and 5, respectively. The bifurcation parameter is λ.
solutions converges to ∞ as λ → 0, see for example [32]. It is therefore not surprising that to this day, no complete rigorous description of the set of equilibria of (2) on Ω = (0, 1) exists. On higher-dimensional domains even the Cahn-Hilliard equilibrium structure is considerably more complicated, see for example the rigorous mathematical results in [18,25], and a complete description seems to be outside the reach of classical mathematical methods.
While rigorous mathematical results are hard to come by, semi-rigorous results are certainly within reach. For example, for the one-dimensional case Ω = (0, 1), the paper [22] explored how bifurcation diagrams for the equilibria of (2) change as a function of σ ≥ 0, based on a combination of rigorous bifurcation-theoretic local studies and numerical path-following. In this paper, the numerical continuation software AUTO [12] was used to determine the bifurcation diagrams based on a spectral method. Two sample bifurcation diagrams are shown in Figure 1, where the left image shows a bifurcation diagram for (2) with µ = 0 and σ = 6, while the right image is for the parameter values µ = 0.3 and σ = 6. In both cases, the vertical axis denotes the L 2 (0, 1)-norm of the solutions. Thus, most of the branches actually represent at least the two different solutions which are related through the reflection symmetry u(x) → u(1 − x). The same diagrams are shown in Figure 2, but now with the vertical axis representing the solution energy (1). In both figures, the solution branches are color-coded by the numerically determined Morse index of the underlying solutions, and the colors black, red, blue, green, magenta, and cyan correspond to indices 0, 1, 2, 3, 4, and 5, respectively.  Figure 2. Numerically computed bifurcation diagrams for the diblock copolymer model (2) on the domain Ω = (0, 1) and with σ = 6. The left diagram is for total mass µ = 0, while the right diagram is for µ = 0.3. The solution measure for both diagrams is the energy of the solutions as defined in (1. The solution branches are color-coded by the numerically determined Morse index of the solutions, and the colors black, red, blue, green, magenta, and cyan correspond to indices 0, 1, 2, 3, 4, and 5, respectively. The bifurcation parameter is λ.  Figure 3. Equilibrium solutions for the diblock copolymer model on the domain Ω = (0, 1) for mass µ = 0 and σ = 6. The left diagram shows solutions for λ = 100, while the right diagram is for λ = 150. In each case, the color of the solution curve indicates its stability as in Figure 1, i.e., black curves are stable solutions, red curves have index 1, and blue functions have index 2.
patterns, many of which are more involved than the ones observed in the Cahn-Hilliard model σ = 0.
Despite the fact that these bifurcation diagrams most likely cannot be established using classical rigorous mathematical techniques, there are other possibilities. One approach is based on using computer-assisted proof techniques, which have for example been used successfully to describe parts of the equilibrium set of the Cahn-Hilliard model on two-dimensional square domains [23]. While the method  Table 1 below, while the two solutions on the right are validated in Table 2. In each case, the color of the solution curve indicates its stability as in Figure 1, i.e., black curves are stable solutions, red curves have index 1, blue functions have index 2, and the green solution has index 3.  Figure 5. Equilibrium solutions for the diblock copolymer model on the domain Ω = (0, 1) for mass µ = 0.3 and σ = 6. The left diagram shows solutions for λ = 100, while the right diagram is for λ = 200. In each case, the color of the solution curve indicates its stability as in Figure 1, i.e., black curves are stable solutions, red curves have index 1, and blue functions have index 2. used in [23] was based on the idea of self-consistent bounds developed in [49], further refinements using the notion of radii polynomials can be found in [10,19,20]. Furthermore, even results concerning heteroclinic connections between equilibrium solutions could be obtained through the use of Conley index techniques [24]. See also the recent results in [44], which should be applicable to finite-dimensional projections of dissipative partial differential equations such as (2).
In the present work, we present a computer-assisted approach to verifying the existence of equilibrium solutions of (2) such as the ones shown in Figures 3, 4, and 5. Rather than using the approach developed in [10,19,20,49], we follow the technique outlined in [35,36,37], see also the extensive literature cited therein. The methods developed in the former papers heavily rely on specific nonlinearity estimates, which are closely tied to the specific solution representation, i.e., to the underlying basis in function space. Even though the estimates have been given in a variety of contexts, for more complicated nonlinearities they appear daunting. In contrast, the method described in [37] relies on standard partial differential equations estimates, which do not depend on the specific solution representation and which should be less affected by the specific form of the nonlinearity. While it is reasonable to expect that this might lead to somewhat less stringent bounds in the computer-assisted proofs, we consider the simplified nonlinearity estimate a worthy trade-off in many situations. The main philosophy of our results follows the outline in [37], yet the necessary estimates for the operator norm of a certain inverse Fréchet derivative is obtained using recent results on eigenvalue exclusion intervals [47]. Morever, in contrast to the finite-element approach that is used in [37] for the construction of the approximative solution, we employ a Galerkin spectral method. On the one hand this allows us to directly use the output from our continuation computations using AUTO, and on the other hand this leads to extremely good solution approximations with very small residuals. Our results are among the first considering a fourthorder partial differential equation using the method of [37], and they are somewhat complicated by the fact that the Fréchet derivative of the right-hand side of (2) is not self-adjoint with respect to the standard L 2 -inner product. For the sake of presentation, we present our results only in the one-dimensional situation Ω = (0, 1) and for the standard double-well potential. However, this method can easily be adapted to other situations and higher dimensions, and this will be presented in the context of Cahn-Morral systems as in [11] in a future paper.
The remainder of this paper is organized as follows. In Section 2 we present an abstract framework for establishing both equilibrium solutions and branches of equilibrium solutions of nonlinear operator equations in Hilbert spaces. The case of single solutions is the subject of Section 2.1, and the case of solution branch pieces is addressed in Section 2.2. While most of the hypotheses for the abstract results can fairly easily be verified in a computer-assisted setting, one exception is a norm estimate for the inverse of an infinite-dimensional operator. For this, we present some recent results on eigenvalue exclusions in Section 2.3. After these preparations, Section 3 is devoted to the special case of the diblock copolymer model. In Section 3.1 we present the functional-analytic framework, which is then used in Section 3.2 to verify all but one of the central assumptions of the abstract validation results. The missing norm estimate for the inverse of the Fréchet derivative of the diblock copolymer operator is then addressed in Section 3.3, where we show how the eigenvalue exclusion results can be used in our setting. After explaining some of the implementation details in Section 3.4, we then present sample computer-assisted proofs in Section 4.
2. An abstract equilibrium validation result. In this section we present an abstract result which guarantees the existence of solutions for a parameter-dependent problem of the form where F : R × X → Y is a Fréchet differentiable nonlinear operator between two Hilbert spaces X and Y . The inner products on these Hilbert spaces are denoted by (·, ·) X and (·, ·) Y , respectively. Our presentation is based on Plum [37], yet slightly modified and extended to fit our applications. In particular, in contrast to the results in [37] we consider the parameter-dependent setting.
Our goal is to establish the existence of a small branch of solutions in a neighborhood of a pair (λ * , u * ) ∈ R × X, where we explicitly do not assume that (λ * , u * ) solves the nonlinear problem (3). Instead, we suppose that the following hypotheses hold.
(H1) The residual of the nonlinear operator F at the pair (λ * , u * ) is small, i.e., there exists a small constant > 0 such that (H2) The Fréchet derivative D u F (λ * , u * ) ∈ L(X, Y ), where L(X, Y ) denotes the Banach space of all bounded linear operators from X into Y , is one-to-one and onto, and its inverse D u F (λ * , u * ) −1 : Y → X is bounded and satisfies Thus, K is an upper bound on the operator norm of D u F (λ * , u * ) −1 ∈ L(Y, X). (H3) For (λ, u) close to (λ * , u * ), the Fréchet derivative D u F (λ, u) ∈ L(X, Y ) satisfies a Lipschitz-type condition. More precisely, there exist positive real constants M 1 , M 2 , d 1 , and d 2 ≥ 0 such that for all pairs (λ, u) ∈ R × X which satisfy both u − u * X ≤ d 1 and |λ − λ * | ≤ d 2 we have where · L(X,Y ) denotes the operator norm in L(X, Y ). (H4) For λ close to λ * , the Fréchet derivative D λ F (λ, u * ) is bounded. More precisely, there exist positive real constants M 3 and d 3 such that for all λ ∈ R with |λ − λ * | ≤ d 3 we have where as usual we identify Y with L(R, Y ).
We would like to point out that all of the constants , K, M k , and d k depend on the choice of the pair (λ * , u * ). Moreover, in our specific application to the diblock copolymer model, we need explicit values for these constants which are determined using rigorous interval computations.

2.
1. An existence result for individual solutions. As we mentioned in the introduction, our goal is to establish the existence of a small solution branch for (3) near the point (λ * , u * ) ∈ R × X. As a first step, we identify conditions on the above four constants which imply that at the fixed parameter value λ * , the nonlinear equation F (λ * , u) = 0 has a unique solution close to u * . This result is a slight reformulation of [37, Theorem 1], but in order to keep the present paper self-contained we include the proof. For this result, we only need hypotheses (H1) and (H2), as well as (H3) for d 2 = 0.

THOMAS WANNER
Then there exists a uniquely determined solution u(λ * ) ∈ X of F (λ * , u) = 0 which satisfies Moreover, this solution u(λ * ) is the unique solution of F (λ * , u) = 0 with In other words, the open X-ball centered at u * and with radius δ 2 contains a unique solution of the nonlinear problem F (λ * , u) = 0, and this solution lies in fact in the closed ball of radius δ 1 centered at u * . Thus, the value of δ 1 quantifies how well the solution u(λ * ) is approximated by the function u * , and the distance δ 2 indicates how isolated the solution u(λ * ) is.
Proof. Consider the quadratic polynomial p(t) = t 2 − 2t/(M 1 K) + 2 /M 1 , whose global minimum is attained at due to assumption (8). The two real roots of p are given by and the second inequality in (8) immediately implies δ 1 < d 1 , as stated in (9). Now let δ be any constant which satisfies δ 1 ≤ δ < δ 2 . Such constants exist since both of the inequalities δ 1 < d 1 and δ 1 <δ 2 hold, and therefore we have δ 1 < δ 2 = min{δ 2 , d 1 }. Furthermore, according to δ < δ 2 ≤δ 2 < δ 3 , the constant δ lies between the two zeros of p, which in turn implies In order to verify the theorem we only need to show that the nonlinear problem F (λ * , u) = 0 has a unique solution in the closed ball Let L = D u F (λ * , u * ) and define the Newton-like nonlinear operator T : X → X by Note that the inverse L −1 exists and is bounded according to (H2). It can easily be seen that the operator T has a fixed point u if and only if u satisfies F (λ * , u) = 0. The remainder of the proof is used to show that T (U δ ) ⊂ U δ , and that T is a uniform contraction -which immediately implies the theorem in combination with Banach's fixed point theorem. For this, let u, v ∈ U δ be arbitrary. Then we have and together with the definition of T , as well as (H2) and (H3) this implies where we used the fact that with u and v also the convex combination , assumption (H3) can be applied. Now let u ∈ U δ be arbitrary. Then the identity T (u * ) = −L −1 F (λ * , u * ) + u * , in combination with the estimate (13) and both (H1) and (H2), yields where for the last inequality we used (11). Thus, T (u) ∈ U δ , as well as T (U δ ) ⊂ U δ .
In order to show that T : U δ → U δ is in fact a contraction, we consider again the estimate in (13), this time for u, v ∈ U δ . Since for all t ∈ [0, 1] we then have tu + (1 − t)v − u * X ≤ δ due to the convexity of U δ , this finally gives According to our choice of δ we have M 1 Kδ < M 1 Kδ 2 ≤ 1, and the proof of the theorem is complete.
Theorem 2.1 establishes the existence of solutions of the nonlinear problem F (λ * , u) = 0 close to an approximate solution u * , and it is similar in spirit to the Newton-Kantorovich theorem, see for example [48,Theorem 5.A]. Its applicability depends of course on the four constants , K, M 1 , and d 1 , which in turn depend on the choice of the pair (λ * , u * ). (Note that M 2 and d 2 have no impact on the result. They will appear in the next theorem.) The constant M 1 can usually be obtained by abstract functional analytic estimates, which normally result in an expression involving certain norms of the approximate solution u * , and the cutoff constant d 1 can frequently be chosen as a fixed number. In fact, in this aspect our result differs slightly from the one in [37]. He considers an estimate which is valid on the entire function space, and which leads to a u-dependent M 1 with higher-order growth as the solution u is far away from u * . In contrast, in our formulation we only allow deviations of u from u * up to distance d 1 , and therefore we can obtain a fixed growth bound M 1 , which in turn allows for explicit formulas for δ 1 and δ 2 . For our application to the diblock copolymer model, we will usually simply choose the constant d 1 = 1. Thus, the constants M 1 and d 1 are generally not that critical.
The situation is different for the remaining constants and K. The former constant is a measure for the accuracy of our numerical approximation, and depending on the numerical method that is used for the construction of u * one should aim to generate as small a residual as possible. We can assume that this can be accomplished if a reasonable numerical method has been chosen. On the other hand, the constant K is a bound on the norm of the inverse operator D u F (λ * , u * ) −1 , and this bound can be large. For example, this will always be the case if (λ * , u * ) is close to a bifurcation point of (3), since the Fréchet derivative D u F necessarily has to be noninvertible at the bifurcation point. If K is large, this needs to be compensated for by a small residual , in order for the estimates in (8) to hold. Notice also that for large values of K the radius δ 2 will be small, i.e., we can only prove uniqueness of the solution u(λ * ) in a small neighborhood of u * . Finally, we would like to stress the fact that among the two estimates in (8), the second one is usually automatically satisfied as long as the first estimate holds.

An existence result for pieces of solution branches.
We now turn our attention to the question of finding pieces of solution branches of the parameterdependent problem (3). Results of this type have already been obtained in [35] in the context of second-order elliptic boundary value problems, combined with a finite element approach. The following abstract result is based on the implicit function theorem and should be seen as a first step towards validating branches. In the future, this has to be further refined by using ideas from branch continuation.
Theorem 2.2 (Branch Existence). Let X and Y be Hilbert spaces, and suppose that the nonlinear parameter-dependent operator F : R×X → Y is Fréchet differentiable. Furthermore, assume that (λ * , u * ) ∈ R × X satisfies hypotheses (H1), (H2), (H3), and (H4). Finally, suppose that In contrast to Theorem 2.1, we now consider the constants δ 1 < δ 2 < δ 3 defined by Then for every choice of In other words, where U δu was defined in (12). The specific value of the constant δ λ depends on the relation between δ u and the constants in (15). More precisely, the constant δ λ can which form the upper boundary of the region of admissible pairs (δ u , δ λ ). The left diagram depicts the normal situation, in which δ 3 < d 1 . In the right image we illustrate the case δ 1 < d 1 < δ 3 , in which the admissible region is reduced due to hypothesis (H3). For both diagrams, we assume that δ λ,opt < d 2 ≤ d 3 , which is usually satisfied.
be chosen as Note that the value δ λ is maximal if we can choose δ u = δ 2 . This is illustrated in Figure 6.
Proof. Due to the first inequality in (14) there exists a constant 0 < q < 1 such that In fact, we can always choose q = 1/2, and the formulation of the theorem uses exactly this choice. However, for later applications we keep the variable q throughout the proof of the result. Let L = D u F (λ * , u * ) ∈ L(X, Y ) denote the Fréchet derivative of F at the approximate solution, which according to (H2) has a continuous inverse, and define the parameter-dependent Newton-like operator Furthermore, for the moment let 0 < δ u ≤ d 1 and 0 < δ λ ≤ min{d 2 , d 3 } be arbitrary constants, and define the closed ball U δu as in (12). Then for all λ ∈ [λ * −δ λ , λ * +δ λ ] and arbitrary u ∈ U δu our assumptions imply In other words, if we now assume that δ u > 0 and δ λ > 0 additionally satisfy Together with the mean value theorem as in [48,Theorem 4.A], this immediately implies be arbitrary. Then due to (H1), (H2), and (H4) we have and in combination with (18) one obtains for all (λ, Thus, the new constraint for the constants δ u > 0 and δ λ > 0 furnishes Together with (18) and 0 < q < 1 this finally implies that T : [λ * − δ λ , λ * + δ λ ] × U δu → U δu is a uniform contraction, and the statement of the theorem follows from the uniform contraction mapping principle, as long as (17) and (19) can be satisfied simultaneously.
For this, we now switch to the special case q = 1/2. One can easily see that equalities in (17) and (19) lead to the blue lines shown in Figure 6, and that the intersections of these lines with the δ u -axis are given by δ 3 and δ 1 , respectively. Moreover, the two lines intersect at δ u = δ 2 , and the corresponding value of δ λ is given by the respective expression in (16). Thus, the light blue region in the left image of Figure 6 shows the maximal admissible set of pairs (δ u , δ λ ). It is nonempty, since the first inequality in (14) is equivalent to δ 1 < δ 3 , the line described by equality in (19) is increasing, and the line associated with (17) is decreasing. Finally, the overall constraints 0 < δ u < d 1 and 0 < δ λ < min{d 2 , d 3 } cannot change this, due to the second estimate in (14). This completes the proof of the theorem.
We would like to point out that the main assumptions in Theorem 2.2 are slightly stronger than the ones in Theorem 2.1, but in return we obtain a piece of a solution branch, rather than just a single solution. Note also that the main condition of the theorem is the first inequality in (14), while the second constraint is mild and usually satisfied.
At first glance, the specific form of the region of admissible (δ u , δ λ )-pairs shown in Figure 6 might seem strange -in particular the fact that the admissible δ λvalues converge to 0 as δ u approaches δ 1 or δ 3 . Notice, however, that for δ u ≈ δ 1 one has to assume that the true solution u(λ * ) will be close to the boundary of the ball U δu , and therefore small perturbations of the parameter λ might lead to the solution u(λ) leaving this fixed ball. Similarly, if we choose δ u ≈ δ 3 , then one has to assume that outside of U δu , but close to its boundary, there might be other nontrivial solutions of (3), and these could enter the ball upon small perturbations of λ.
Finally, the formula for the optimal, i.e., largest, value of δ λ given in (16) could be improved. For this, recall that while in the proof of Theorem 2.2 we introduced the Lipschitz constant 0 < q < 1, at the end we only considered the special case q = 1/2. One can in fact reformulate the statement of the theorem using q directly, and thereby obtain a q-dependent formula for the optimal δ λ -value. This value can then be maximized by choosing an appropriate value of q.

Exclusion of eigenvalues.
One of the central issues in the approach presented so far is the ability to obtain good estimates for the constant K in hypothesis (H2). For this, we make use of recent results due to Watanabe, Nagatou, Plum, and Nakao [47], which build upon original work by Nagatou and collaborators [30,31]. In their paper, two approaches are developed for finding regions in the complex plane in which certain generalized eigenvalue problems have no eigenvalues. While their methods are formulated in a very general setting, we only need a small subset of the results in [47]. Thus, in order to keep the present paper self-contained, we include this reduced version in the remainder of this section. We consider a generalized eigenvalue problem of the form where the linear operators A and B satisfy the following assumptions.
(H5) Let X ⊂ Y denote two Hilbert spaces with inner products (·, ·) X and (·, ·) Y . Suppose that A : D(A) ⊂ X → Y is a (generally unbounded) linear operator which maps its domain D(A) bijectively onto Y, and assume that its (H6) With the spaces X and Y as before, the linear operator B : X → Y is bounded, i.e., there exists a constant M 4 ≥ 0 such that (H7) For every integer N ∈ N there exists an N -dimensional subspace X N ⊂ X such that the following holds. If P N ∈ L(X , X ) denotes the orthogonal projection onto X N , then there exists a constant C N > 0 with lim Hypotheses (H5), (H6), and (H7) form the theoretical backbone for the eigenvalue exclusion result below. However, we need another assumption which is very much numerical at heart. Since we are interested in excluding eigenvalues, we will in one way or another have to establish the invertibility of a certain operator. This has to be done in a finite-dimensional setting, i.e., we will need to guarantee the invertibility of certain matrices, whose building blocks will be presented next. As (H7) indicates, we intend to treat the infinite-dimensional eigenvalue problem (20) through finite-dimensional approximations of the Hilbert space X . Thus, let N ∈ N be arbitrary, let X N be as in (H7), and let {ϕ 1 , . . . , ϕ N } denote a basis of X N . Then we define the N × N -matrices Q, R, B ∈ R N ×N via where as usual m, n = 1, . . . , N denote the row and column indices of the entries, respectively. One can easily see that the matrices Q and R are symmetric and positive definite. Therefore, each of them has a Cholesky factorization, i.e., one has Q =QQ T and R =RR T for lower triangularQ,R ∈ R N ×N . (25) Then our last hypothesis can be formulated as follows.
(H8) For N ∈ N let {ϕ 1 , . . . , ϕ N } denote a basis of the space X N from (H7), and consider the matrices defined in (24) and (25). Furthermore, let γ 0 be an arbitrary, but fixed, real number. Then we assume that the matrix Q − γ 0 B is invertible, and that the positive constant τ is chosen such that where · 2 denotes the standard matrix 2-norm, i.e., its largest singular value. Then the following result holds.

Theorem 2.3 (Eigenvalue Exclusion).
Consider the generalized eigenvalue problem (20), and suppose that γ 0 ∈ R is chosen in such a way that (H5), (H6), (H7), and (H8) are satisfied. Finally, assume that N ∈ N is such that Then (20) has no eigenvalue in the interval (γ 0 − γ 1 , γ 0 + γ 1 ), where Proof. Our verification of this result closely follows the proof of [47, Theorem 5.1]. Thus, we begin by considering arbitrary elements v ∈ D(A) and w ∈ Y which satisfy Let X N be the finite-dimensional subspace of X from (H7), and let X ∞ denote its orthogonal complement in X . Then we can decompose v as and after projection onto X N and X ∞ the last identity in (29) can equivalently be written as In other words, we have the identity (w ∞ , ϕ m ) Y = 0 for all m = 1, . . . , N . Now let m ∈ {1, . . . , N } be arbitrary. Together with (21), (24), and (30), these definitions and the orthogonality properties of w ∞ imply the identities and the definition of the matrix product furnishes Rw = (Q − γ 0 B)v. Now let · 2 denote the standard Euclidean norm in R N . Then (24) and (H8) further imply since γ 0 Bv ∞ +w−w ∞ ∈ X N and w ∞ is Y-orthogonal to X N . An application of (H6) now finally gives We now turn our attention to bounding v ∞ . According to (31), (H6), and (H7) one obtains where for the last inequality we used (32). But now (27) implies and with (32) one further gets Combining the above two estimates furnishes where γ 1 was defined in (28). Recalling our definition of v and w in (29), we therefore have proved that and this immediately implies that the operator A − γ 0 B : D(A) → Y is one-to-one. Consider now the operator F = I − γ 0 A −1 B ∈ L(X , X ). If v ∈ X lies in its nullspace, then we have v = γ 0 A −1 Bv ∈ D(A). But in this case Av − γ 0 Bv = AFv = 0, and the just established injectivity of A − γ 0 B therefore yields v = 0. In other words, F ∈ L(X , X ) is one-to-one. Notice, however, that due to (H5) the operator A −1 B ∈ L(X , X ) is compact, and applying the Fredholm alternative to F then implies that F −1 ∈ L(X , X ) exists. Now let w ∈ Y be arbitrary. Then we have A −1 w ∈ D(A) ⊂ X , and therefore there exists an element v ∈ X with A −1 w = Fv = v − γ 0 A −1 Bv. This implies that v ∈ D(A), as well as the identity Av − γ 0 Bv = w. In other words, the operator A − γ 0 B : D(A) → Y is one-to-one and onto.
After these preparations we are finally in a position to establish the eigenvalue exclusion result. For this, let v ∈ D(A) \ {0} be such that Av = γBv for some real number γ ∈ R, i.e., the number γ is an eigenvalue of (20). Then (33) implies in combination with (H6) the estimate and therefore |γ − γ 0 | ≥ γ 1 . This completes the proof of the theorem.
In our applications later on, this result will be applied multiple times, so that we can cover a large eigenvalue exclusion interval and obtain reasonable estimates for the constant K in (H2).
3. Application to the diblock copolymer model. In this section we demonstrate how the diblock copolymer model (2) can be studied using the theory presented in the previous section. For this, we do not precisely consider the form of the model given in the introduction. Rather, we incorporate the total mass value µ into the equation by considering the transformation u → µ + u. In this way, we have the mass constraint for u fixed at zero. If we further define the nonlinearity f as f (u) = −W (u), which for the standard double-well potential W implies the identity f (u) = u − u 3 , then we consider the nonlinear operator equation Note that even though both µ and σ are parameters of the model, we only emphasize the dependence on λ in the definition of F , as this is our main bifurcation parameter.
3.1. Functional-analytic framework. In the following, we concentrate on functions u : Ω → R defined on the one-dimensional domain Ω = (0, 1). Our primary function spaces will be related to the standard Sobolev spaces H k (0, 1) = W k,2 (0, 1), see for example [1]. Particularly important is the space H 1 (0, 1) of all weakly differentiable L 2 (0, 1)-functions whose weak derivative u is square-integrable, equipped with the standard norm Due to the zero mass constraint, we specifically consider subspaces H k of the Sobolev spaces H k (0, 1), which are defined as follows. Note that the eigenvalues and eigenfunctions of the negative Laplacian −∆ on the one-dimensional domain Ω = (0, 1) and subject to homogeneous Neumann boundary conditions, are given by the constant function ϕ 0 (x) = 1 with eigenvalue κ 0 = 0, together with the sequence κ = 2 π 2 and ϕ (x) = √ 2 cos πx for ∈ N .
Then we define Note that in the definition of H k , the constant Laplacian eigenfunction ϕ 0 is omitted due to the zero mass constraint. For k ≥ 0 the coefficients α are given by α = (u, ϕ ) L 2 (0,1) , since the eigenfunctions ϕ form a complete orthonormal set in L 2 (0, 1). For k > 0 and u ∈ H −k the series is interpreted formally, i.e., the element u ∈ H −k is identified with the sequence of its Fourier coefficients, and one can easily see that in this case u acts as a bounded linear functional on H k . In fact, for every k ∈ Z the space H k is a closed subspace of the Sobolev space H k (0, 1). Furthermore, H k is a Hilbert space with respect to the inner product (36) Notice that for example the space H 2 consists of all functions in H 2 (0, 1) which satisfy homogeneous Neumann boundary conditions and whose integral over the domain Ω = (0, 1) vanishes. Moreover, the product (u, v) H k is also defined for functions u ∈ H p and v ∈ H q as long as 2k ≤ p + q, since The above choice of spaces and norms has some convenient side-effects. First of all, one can easily see that the Laplacian ∆ is in fact an isometry between the spaces H k and H k−2 for all k ∈ Z. In other words, we have In this identity we slightly abuse notation by abbreviating (∆| H k+2 ) −1 by ∆ −1 . Moreover, we can define the square root (−∆) 1/2 of the negative Laplacian via which is an isometry between H k and H k−1 and satisfies In addition, for all m ∈ Z the identities and are satisfied, as long as all inner products are defined. Finally, we can obtain a number of explicit norm bounds which are necessary for the implementation of the method described in the previous section. These bounds are collected in the following lemma.
(a) If · ∞ denotes the maximum norm, then for all u ∈ H 1 (0, 1) ⊂ C[0, 1] we have 1) are arbitrary, then also u · v ∈ H 1 (0, 1) and we have H 1 (0, 1) is arbitrary, then we have (d) If u ∈ H k is arbitrary, and if m < k are arbitrary integers, then u ∈ H m and we have Proof. The estimate in (a) has been established in [28,38,46], where it is also shown that the constant C a is optimal, with equality holding for u(x) = cosh x. For the inequality in (b) one obtains where we use the estimate established in (a) and abbreviate L 2 (0, 1) to Then one can easily see that and the inequality u H 1 ≤ u H 1 (0,1) follows similarly. Finally, if u has the expansion above and is contained in H k , then which establishes (d).
The above estimates are central for our discussions below. In addition, we will need to be able to project elements in the Sobolev space H 1 (0, 1) into H 1 . This will be accomplished by the projection P : Note that every function w ∈ H 1 (0, 1) has a series representation w = ∞ =0 α ϕ where α = (w, ϕ ) L 2 (0,1) . Then one can easily see that P w = ∞ =1 α ϕ . Furthermore, we have

THOMAS WANNER
3.2. Fundamental nonlinearity estimates. After these preparations, we finally return to the diblock copolymer equation (2). For this, let µ ∈ R and σ ∈ R + 0 be arbitrary, but fixed, real constants. Then we consider the nonlinear operator F defined by (42) and where for the sake of presentation we always assume that f (u) = u−u 3 , even though this can easily be changed. It follows from standard results that this operator is Fréchet differentiable, and the following proposition demonstrates that (H3) is satisfied.
Proposition 3.2 (Verification of (H3)). Let µ ∈ R and σ ∈ R + 0 be constants, and let f (u) = u − u 3 . Furthermore, consider the nonlinear operator F defined in (42), acting on the spaces X = H 1 and Y = H −3 defined in (35). Finally, let d 1 > 0 and d 2 ≥ 0 be arbitrary constants and let (λ * , u * ) ∈ R × X. Then for all , and In these formulas, · L(X,Y ) denotes the operator norm in L(X, Y ), and the constants C a and C c are as in Lemma 3.1. Note that d 1 > 0 and d 2 ≥ 0 can be chosen arbitrarily, but their specific choice affects M 1 .
Proof. Consider the projection P : H 1 (0, 1) → H 1 defined in (40). Then for all functions v ∈ X = H 1 and (λ,ū) ∈ R × X the Fréchet derivative satisfies For all v ∈ X with v X ≤ 1 this implies, together with Lemma 3.1 and (37), the estimate where we also used P w H 1 ≤ P w H 1 (0,1) ≤ w H 1 (0,1) for all w ∈ H 1 (0, 1). The above estimate readily implies We continue by first bounding the difference λf (µ+u)−λf (µ+u * ) in the H 1 (0, 1)norm. According to f (u) = 1 − 3u 2 and Lemma 3.1 one obtains and C b = 2C a , one can easily see that F satisfies (H3) with the constants M 1 and M 2 as in the formulation of the proposition.
The verification of hypothesis (H4) is the subject of the next result.
Proposition 3.3 (Verification of (H4)). Let µ ∈ R and σ ∈ R + 0 be constants, and let f (u) = u − u 3 . Furthermore, consider the nonlinear operator F defined in (42), acting on the spaces X = H 1 and Y = H −3 defined in (35). Then for all λ ∈ R and u * ∈ X we have In this formula, we identify Y with L(R, Y ), and P : H 1 (0, 1) → H 1 denotes the projection defined in (40). In other words, the nonlinear operator F satisfies hypothesis (H4) with the above constant M 3 and with d 3 = ∞.

3.3.
Bounding the norm of the inverse Fréchet Derivative. We now turn our attention to hypothesis (H2), i.e., the bound on the operator norm for the inverse Fréchet derivative. For this, let µ ∈ R and σ ∈ R + 0 again be arbitrary, but fixed, real constants. In addition, let λ * ∈ R and u * ∈ H 1 be arbitrary. Then the nonlinear operator F defined in (42) is Fréchet differentiable with respect to u at the point (λ * , u * ), and if we denote this derivative by L = D u F (λ * , u * ), then we have and where as before we only consider f (u) = u − u 3 . Then the following result shows that one can obtain an explicit bound on the operator norm of L −1 by solving an associated self-adjoint eigenvalue problem. Note that the eigenvalue problem below differs from the one considered in [3]. Moreover, our approach differs from the one described in [36,Section 3.3], which is based on the method by Lehmann and Goerisch and requires a homotopy argument. The latter is not necessary in our approach, as we rely on eigenvalue exclusions results.

Proposition 3.4 (Norm Estimate via Eigenvalue Exclusion).
Consider µ ∈ R and choose σ ∈ R + 0 , define f (u) = u − u 3 , and let λ * ∈ R and u * ∈ H 1 be arbitrary. Furthermore, consider the bounded linear operator L defined in (43), acting between the spaces X = H 1 and Y = H −3 defined in (35). Finally, consider the bounded linear operator S ∈ L(X, X) defined by Sv = (−∆) −2 Lv , (44) and assume that there exists a constant K > 1 such that S has no eigenvalue in the interval Then the operator L has a bounded inverse L −1 ∈ L(Y, X) and we have In other words, hypothesis (H2) is satisfied.
Proof. The properties of L and the fact that ∆ : H k → H k−2 is an isometry immediately imply that S ∈ L(X, X). Furthermore, notice that S can be rewritten in the form where the projection P : H 1 (0, 1) → H 1 was defined in (40). Note that for every function v ∈ H 1 the second and third terms in this representation of S satisfy Due to the compact embeddings of both H 3 and H 5 in the space X = H 1 we therefore see that S is a compact perturbation of the negative identity in X.
As our second step, we now verify that the operator S ∈ L(X, X) is symmetric, and therefore self-adjoint. For this, let v, w ∈ X be arbitrary. Then (38), (39), since the second to last expression is completely symmetric in v and w. So far we have shown that S ∈ L(X, X) is a self-adjoint compact perturbation of the negative identity on the separable Hilbert space X = H 1 . Then [9, Theorem 4.2.23] guarantees the existence of a complete orthonormal set {ψ } ∞ =1 in X consisting of eigenfunctions of S, i.e., there exist real numbers ν which converge to −1 as → ∞ and such that Sψ = ν ψ for all ∈ N. Moreover, each eigenvalue ν = −1 has finite multiplicity. Now let w ∈ X be arbitrary, and define α = (w, ψ ) X for ∈ N. Then According to (45) we have |ν | ≥ 1/K for all ∈ N, and due to the X-orthonormality of the ψ one obtains In combination with (37) this implies for all w ∈ X, and since L is a Fredholm operator with index zero we therefore have the existence of L −1 ∈ L(Y, X). The estimate (46) finally follows if we let w = L −1 v for v ∈ Y . This completes the proof of the proposition.
The above proposition shows that obtaining an explicit bound on the operator norm of the inverse Fréchet derivative can be reduced to an eigenvalue exclusion problem. At first glance, the choice of the operator S seems somewhat arbitrary. However, this choice is exactly necessary if we want to obtain precise estimates for the operator norm of L −1 . To see this, let 0 = w ∈ Y = H −3 and v ∈ X = H 1 be related via the identity Lv = w. Then we have If on the left hand side we take the supremum over all w ∈ Y \ {0}, then we obtain exactly the operator norm of L −1 in L(Y, X). Moreover, this supremum equals the supremum of the expression on the right over v ∈ X \ {0}, which in turn is the same as the reciprocal of the infimum of Sv X / v X . According to the proof of Proposition 3.4, this implies In other words, the eigenvalue problem for S determines the operator norm of L −1 precisely. Yet, since the operator S defined in (44) is not in a form which is directly amenable to the results of Section 2.3, we need to introduce a slight reformulation. This is the subject of the next result.
Proposition 3.5 (Verification of (H2)). Let µ ∈ R and σ ∈ R + 0 , let f (u) = u − u 3 , and let λ * ∈ R and u * ∈ H 1 be arbitrary. Furthermore, consider the generalized eigenvalue problem where we define Finally, let K > 1 be a constant such that (48) has no eigenvalue γ in the interval Then the linear operator L defined in (43) has a bounded inverse L −1 ∈ L(Y, X) and we have In other words, hypothesis (H2) is satisfied.
Proof. Consider the operator S defined in (44), and suppose that η ∈ R \ {−1} is an eigenvalue for S with associated eigenfunction v ∈ H 1 \ {0}. Then we have From the last identity one obtains v = −A −1 Bv/(1 + η) ∈ H 3 , since B : H 1 → H −1 and A is an isometry between H 3 and H −1 .
The setting of Proposition 3.5 enables us to use the result on eigenvalue exclusions, i.e., the use of Theorem 2.3. For this, hypotheses (H5) through (H8) have to be verified, and the following result addresses the first two of these. Proposition 3.6 (Verification of (H5) and (H6)). Let µ ∈ R and σ ∈ R + 0 be given, consider the nonlinearity f (u) = u − u 3 , and let λ * ∈ R and u * ∈ H 1 be arbitrary. Furthermore, consider the generalized eigenvalue problem Av = γBv, where Av = (−∆) 2 v and Bv = λ * ∆(f (µ + u * )v) + λ * σv. If we consider the spaces then hypothesis (H5) is satisfied. Moreover, for all v ∈ X one has i.e., hypothesis (H6) holds as well.
Proof. The inclusions D(A) ⊂ X ⊂ Y are obvious, and since the embedding D(A) → X is compact, the compactness of A −1 ∈ L(Y, X ) follows from the fact that A −1 : H −1 → H 3 is an isometry. Now let v ∈ H 3 and w ∈ H 1 be arbitrary, and similarly to (36) In other words, (H5) is satisfied. As for hypothesis (H6), consider the projection P : H 1 (0, 1) → H 1 defined in (40). Then for arbitrary v ∈ H 1 Lemma 3.1 and (37) where we again used P w H 1 ≤ P w H 1 (0,1) ≤ w H 1 (0,1) for all w ∈ H 1 (0, 1). This immediately implies (H6) with the constant M 4 as in the formulation of the proposition.
For the next hypothesis we need to define appropriate finite-dimensional subspaces X N ⊂ X . As our approach is based on spectral methods, we make use of the functions ϕ defined in (34). This leads to the following result.
Proposition 3.7 (Verification of (H7)). Consider the functions ϕ defined in (34), which form a complete orthonormal set in H 0 , and let N ∈ N. Define the Ndimensional subspace Furthermore, let A = (−∆) 2 be as in Proposition 3.6, and let P N ∈ L(X , X ) denote the orthogonal projection onto X N . Then for all v ∈ D(A) we have In other words, (H7) holds.
This completes the proof. At this point, we only have to worry about hypothesis (H8). As we mentioned earlier, this hypothesis is essentially numerical at heart, and we could leave its discussion to the next section. However, for our specific choice of the spaces X N , it turns out that the matrices Q and R are diagonal matrices, and therefore we can simplify the formulation of (H8). This is the subject of the following final result.
Proof. We begin by determining the matrices Q and R defined in (24). If δ m,n denotes the Kronecker delta, then X = H 1 and Y = H −1 imply for m, n = 1, . . . , N the identities Q m,n = (ϕ n , ϕ m ) X = δ m,n κ m and R m,n = (ϕ n , ϕ m ) Y = δ m,n κ −1 m , i.e., both Q and R are diagonal matrices. This furnisheŝ As for the matrix B defined in (24), the definition of B in Proposition 3.6, together with identities (38), (39), (40), and (41), yields Now define S =QQR and T =QBR. Then for m, n = 1, . . . , N we have the identities as well as which is exactly how these matrices were defined in (51). Now assume that the matrix S − γ 0 T is invertible and that (52) holds. According to the definitions of S and T we have From this, the result follows immediately.
3.4. Implementation details. The propositions of the last two subsections presented explicit formulas for all constants appearing in our eight main hypotheses. These formulas depend on the approximate solution pair (λ * , u * ) ∈ R × H 1 , and in the following we always assume that the function u * is given explicitly via the sum where ϕ was defined in (34). In fact, the numerical bifurcation diagrams shown in the introduction were computed using a spectral method, and therefore provide exactly the coefficients α in this sum. The specific choice of M depends both on the solution branch and the parameter range that is being considered. If the approximate solution u * is given as in (53), and since we only consider the polynomial nonlinearity f (u) = u − u 3 , the necessary H k -norms of both f (µ + u * ) and f (µ + u * ) can now easily be determined from the Fourier coefficients α , using standard convolution product techniques. In order to use these computations in a computer-assisted proof setting, one of course has to take roundoff errors into account. This is accomplished by using interval arithmetic, as described for example in [29]. In our demonstrations below, we make use of the Matlab toolbox Intlab [39].
While the norm computations mentioned in the last paragraph can easily be performed in interval arithmetic, the verification of hypothesis (H8) requires more work. Suppose we are in the situation of Proposition 3.8. Then the matrix S − γ 0 T will in fact be given as an interval matrix, i.e., every matrix entry is a small interval of numbers, and S − γ 0 T is the set Ξ of all matrices whose entries lie in these intervals. In order to establish a computer-assisted proof, one has to show the invertibility of all matrices in Ξ, and similarly, (52) has to hold for all of them. This can be accomplished using the following technique from [41]. Let J ∈ R N ×N be the numerically computed inverse of an arbitrary matrix in Ξ, for example the matrix whose entries are exactly the centers of the entry intervals. Then we can use interval computations to find rigorous bounds 0 < 1 < 1 and 2 > 0 with Then [41, Lemma 2.1] implies that every matrix in Ξ is invertible with In other words, the only computations that have to be performed in interval arithmetic are matrix multiplications and additions, as well as finding tight bounds on the matrix 2-norm of an interval matrix -and for the latter Intlab provides sophisticated algorithms. For more details we refer the reader to the survey [40].
4. Computational results. In this final section we combine everything to present sample computer-assisted proofs for the diblock copolymer model. We mostly concentrate on the validation of individual solutions via Theorem 2.1, and only briefly illustrate the second Theorem 2.2 which concentrates on perturbation robustness. All computations consider the diblock copolymer model on the one-dimensional domain Ω = (0, 1) with parameter σ = 6, and both for mass µ = 0 and µ = 0.3. More specifically, our goal is to rigorously establish the existence of some of the equilibrium solutions shown in the bifurcation diagrams in Figures 1 and 2. They were performed using Matlab and Intlab [39].
Equilibrium Validation for µ = 0 with K = 500 λ Number Index δ 1 δ 2 N 200 # 2 0 2.642564e-14 1.321300e-11 4.858832e-07 1828 # 3 1 6.084200e-14 3.042193e-11 4.959956e-07 1108 Table 2. Rigorous numerical results for the diblock copolymer model with µ = 0. For λ = 200 we consider representative solutions from the second and third branches shown in the left diagram of Figure 1, numbered from top to bottom. The table lists the index of each solution, the value in (H1), as well as δ 1 and δ 2 from Theorem 2.1. The last column indicates which value of N is chosen in Theorem 2.3. The computations use K = 500.
We first turn our attention to the case µ = 0, i.e., to the left diagram of Figure 1. For each of the three λ-values 100, 150, and 200 we consider representative solutions from each of the branches shown in this numerically determined bifurcation diagram. We number these solutions according to their L 2 -norm, i.e., from top to bottom in the bifurcation diagram. Since the bifurcation diagram was computed by AUTO [12] using a spectral method based on (53) with M = 100, we used these solution representations as approximations for u * . We then tried to validate the solutions using Theorem 2.1 by picking K = 100 and K = 500 in (H2) and Proposition 3.5, and the results of these computations are contained in Tables 1 and 2, respectively. We would like to point out that due to the spectral representation of the approximate solutions u * we are able to achieve extremely small values of the residual . With K = 100, all but two solutions could be rigorously validated with L ∞ -norm errors of size 10 −12 , see Table 1. The two exceptional cases, i.e., the solutions on the second and third branch from the top for λ = 200, took a value of K = 500 to validate, as seen in Table 2. These solutions are shown in the right image of Figure 4, and their shape is similar, i.e., they are fairly close to each other in function space. This results in higher norms of the inverse of the Fréchet derivatives. In these cases, the obtained error bounds are of order 10 −11 . In other words, these computations prove that up to extremely small errors, the diblock copolymer model for µ = 0 has equilibria as the ones shown in Figures 3 and 4 in the introduction.
Similar results can be obtained for µ = 0.3, see Table 3. In this case we consider the λ-values 100 and 200, and again representatives from each of the branches shown in the right diagram of Figure 1 numbered from top to bottom, i.e., according to decreasing L 2 -norm. Notice that now the value K = 100 works for all nine equilibrium solutions. Again, these computations prove that up to small errors, the diblock copolymer model for µ = 0.3 has equilibria such as the ones shown in Figure 5 in the introduction.
of the branches shown in the left diagram of Figure 1, numbered from top to bottom. The table lists the values of δ 1 , δ 2 , δ 3 , and δ λ,opt from Theorem 2.2. These values show that in all cases we were able to prove the existence of a small solution branch close to the approximate solution pair (λ * , u * ). The last two columns indicate the numerical estimate K est for the optimal constant K in (H2), as well as the value of K used for the validation computation.