Computer-assisted proof of heteroclinic connections in the   one-dimensional Ohta-Kawasaki model

Jacek Cyranka; Thomas Wanner

arXiv:1703.01022·math.AP·May 29, 2020

Computer-assisted proof of heteroclinic connections in the one-dimensional Ohta-Kawasaki model

Jacek Cyranka, Thomas Wanner

PDF

1 Repo

TL;DR

This paper uses computer-assisted methods to rigorously prove the existence of heteroclinic connections in a complex PDE model of diblock copolymers, revealing intricate dynamical behaviors not seen in simpler related models.

Contribution

The paper provides the first rigorous proof of heteroclinic connections in the one-dimensional Ohta-Kawasaki model, highlighting its multistability and complex attractor structure.

Findings

01

Existence of heteroclinic connections between homogeneous and energy-minimizing states.

02

Not all solutions near the homogeneous state converge to the global minimizer.

03

Presence of stable states with higher energy trapping solutions.

Abstract

We present a computer-assisted proof of heteroclinic connections in the one-dimensional Ohta-Kawasaki model of diblock copolymers. The model is a fourth-order parabolic partial differential equation subject to homogeneous Neumann boundary conditions, which contains as a special case the celebrated Cahn-Hilliard equation. While the attractor structure of the latter model is completely understood for one-dimensional domains, the diblock copolymer extension exhibits considerably richer long-term dynamical behavior, which includes a high level of multistability. In this paper, we establish the existence of certain heteroclinic connections between the homogeneous equilibrium state, which represents a perfect copolymer mixture, and all local and global energy minimizers. In this way, we show that not every solution originating near the homogeneous state will converge to the global energy…

Equations139

E_{ϵ, σ} [u] = \int_{Ω} (\frac{ϵ ^{2}}{2} ∣\nabla u ∣^{2} + F (u)) d x + \frac{σ}{2} \int_{Ω} (- Δ)^{- 1/2} (u (x) - μ)^{2} d x .

E_{ϵ, σ} [u] = \int_{Ω} (\frac{ϵ ^{2}}{2} ∣\nabla u ∣^{2} + F (u)) d x + \frac{σ}{2} \int_{Ω} (- Δ)^{- 1/2} (u (x) - μ)^{2} d x .

u_{t}

u_{t}

u_{t} = - (u_{xx} + λ f (u))_{xx} - λσ (u - μ),

u_{t} = - (u_{xx} + λ f (u))_{xx} - λσ (u - μ),

μ = \frac{1}{L} \int_{0}^{L} u (t, x) d x \mbox an d u_{x} (t, x) = u_{xxx} (t, x) = 0 \mbox f or t > 0, x \in {0, L} .

μ = \frac{1}{L} \int_{0}^{L} u (t, x) d x \mbox an d u_{x} (t, x) = u_{xxx} (t, x) = 0 \mbox f or t > 0, x \in {0, L} .

u (t, x) = a_{0} (t) + 2 k = 1 \sum \infty a_{k} (t) cos \frac{π k x}{L},

u (t, x) = a_{0} (t) + 2 k = 1 \sum \infty a_{k} (t) cos \frac{π k x}{L},

a_{0} (t) = 0 \mbox f or a l l t \geq 0 .

a_{0} (t) = 0 \mbox f or a l l t \geq 0 .

e_{k} (x) := cos \frac{π k x}{L}, \mbox w h er e k \geq 0,

e_{k} (x) := cos \frac{π k x}{L}, \mbox w h er e k \geq 0,

u^{3}

u^{3}

u^{3} = j_{1}, j_{2}, j_{3} \in Z j_{1} + j_{2} + j_{3} = 0 \sum a_{j_{1}} a_{j_{2}} a_{j_{3}} + 2 k \in N \sum j_{1}, j_{2}, j_{3} \in Z j_{1} + j_{2} + j_{3} = k \sum a_{j_{1}} a_{j_{2}} a_{j_{3}} cos \frac{π k x}{L}

u^{3} = j_{1}, j_{2}, j_{3} \in Z j_{1} + j_{2} + j_{3} = 0 \sum a_{j_{1}} a_{j_{2}} a_{j_{3}} + 2 k \in N \sum j_{1}, j_{2}, j_{3} \in Z j_{1} + j_{2} + j_{3} = k \sum a_{j_{1}} a_{j_{2}} a_{j_{3}} cos \frac{π k x}{L}

\frac{d a _{k}}{d t} = (- \frac{k ^{4} π ^{4}}{L ^{4}} + \frac{λ k ^{2} π ^{2}}{L ^{2}} - λσ) a_{k} - \frac{λ k ^{2} π ^{2}}{L ^{2}} j_{1}, j_{2}, j_{3} \in Z j_{1} + j_{2} + j_{3} = k \sum a_{j_{1}} a_{j_{2}} a_{j_{3}} for all k > 0 .

\frac{d a _{k}}{d t} = (- \frac{k ^{4} π ^{4}}{L ^{4}} + \frac{λ k ^{2} π ^{2}}{L ^{2}} - λσ) a_{k} - \frac{λ k ^{2} π ^{2}}{L ^{2}} j_{1}, j_{2}, j_{3} \in Z j_{1} + j_{2} + j_{3} = k \sum a_{j_{1}} a_{j_{2}} a_{j_{3}} for all k > 0 .

a_{0} = 0, \mbox a s w e l l a s a_{k} = a_{- k} \in R \mbox f or a l l k \in Z .

a_{0} = 0, \mbox a s w e l l a s a_{k} = a_{- k} \in R \mbox f or a l l k \in Z .

H := {{a_{k}}_{k \in Z} \in H : there exists a constant C \geq 0 such that ∣ a_{k} ∣ \leq \frac{C}{∤ k ∤ ^{6}} for k \in Z},

H := {{a_{k}}_{k \in Z} \in H : there exists a constant C \geq 0 such that ∣ a_{k} ∣ \leq \frac{C}{∤ k ∤ ^{6}} for k \in Z},

H^{'} := H \cap {{a_{k}}_{k \in Z} : a_{0} = 0 and a_{k} = a_{- k} for all k \in Z} .

H^{'} := H \cap {{a_{k}}_{k \in Z} : a_{0} = 0 and a_{k} = a_{- k} for all k \in Z} .

\frac{d a _{k}}{d t} = (- k^{2} (\frac{π}{L})^{2} (k^{2} (\frac{π}{L})^{2} - λ) - λσ) a_{k} - λ k^{2} (\frac{π}{L})^{2} ∣ j_{1} ∣, ∣ j_{2} ∣, ∣ j_{3} ∣ \leq n j_{1} + j_{2} + j_{3} = k \sum a_{j_{1}} a_{j_{2}} a_{j_{3}},

\frac{d a _{k}}{d t} = (- k^{2} (\frac{π}{L})^{2} (k^{2} (\frac{π}{L})^{2} - λ) - λσ) a_{k} - λ k^{2} (\frac{π}{L})^{2} ∣ j_{1} ∣, ∣ j_{2} ∣, ∣ j_{3} ∣ \leq n j_{1} + j_{2} + j_{3} = k \sum a_{j_{1}} a_{j_{2}} a_{j_{3}},

[x] = k = 1 \prod n [x_{k}^{-}, x_{k}^{+}] \subset R^{n}, \mbox a s l o n g a s x_{k}^{-} \leq x_{k}^{+} \mbox f or k = 1, \dots, n .

[x] = k = 1 \prod n [x_{k}^{-}, x_{k}^{+}] \subset R^{n}, \mbox a s l o n g a s x_{k}^{-} \leq x_{k}^{+} \mbox f or k = 1, \dots, n .

[df (Z)] = {M = (M_{ij})_{i, j = 1}^{n} \in R^{n \times n} M_{ij} \in [z \in Z in f \frac{\partial f _{i}}{\partial x _{j}} (z), z \in Z sup \frac{\partial f _{i}}{\partial x _{j}} (z)]} .

[df (Z)] = {M = (M_{ij})_{i, j = 1}^{n} \in R^{n \times n} M_{ij} \in [z \in Z in f \frac{\partial f _{i}}{\partial x _{j}} (z), z \in Z sup \frac{\partial f _{i}}{\partial x _{j}} (z)]} .

z^{'} = f (z), z \in R^{n}, f \in C^{2} (R^{n}, R^{n}) .

z^{'} = f (z), z \in R^{n}, f \in C^{2} (R^{n}, R^{n}) .

c_{N} (N) = \overline{B_{u (N)}} \times \overline{B_{s (N)}} .

c_{N} (N) = \overline{B_{u (N)}} \times \overline{B_{s (N)}} .

\begin{array}[]{lclclclclcl}\displaystyle N_{c}^{-}&:=&\displaystyle\partial{B_{u(N)}}\times\overline{B_{s(N)}}\;,&&\displaystyle N_{c}^{+}&:=&\displaystyle\overline{B_{u(N)}}\times\partial B_{s(N)}\;,&&\displaystyle N_{c}&:=&\displaystyle\overline{B_{u(N)}}\times\overline{B_{s(N)}}\;,\\[8.61108pt] \displaystyle N^{-}&:=&\displaystyle c_{N}^{-1}\left(N_{c}^{-}\right)\;,&&\displaystyle N^{+}&:=&\displaystyle c_{N}^{-1}\left(N_{c}^{+}\right)\;,&&N&:=&c_{N}^{-1}\left(N_{c}\right).\end{array}

\begin{array}[]{lclclclclcl}\displaystyle N_{c}^{-}&:=&\displaystyle\partial{B_{u(N)}}\times\overline{B_{s(N)}}\;,&&\displaystyle N_{c}^{+}&:=&\displaystyle\overline{B_{u(N)}}\times\partial B_{s(N)}\;,&&\displaystyle N_{c}&:=&\displaystyle\overline{B_{u(N)}}\times\overline{B_{s(N)}}\;,\\[8.61108pt] \displaystyle N^{-}&:=&\displaystyle c_{N}^{-1}\left(N_{c}^{-}\right)\;,&&\displaystyle N^{+}&:=&\displaystyle c_{N}^{-1}\left(N_{c}^{+}\right)\;,&&N&:=&c_{N}^{-1}\left(N_{c}\right).\end{array}

h (0, x) = b_{c} (x)

h (0, x) = b_{c} (x)

h (1, x) = (x, 0)

h (t, x) \in N_{c}^{-}

Q (x, y) = α (x) - β (y), \mbox f or a r bi t r a r y (x, y) \in R^{u (N)} \times R^{s (N)},

Q (x, y) = α (x) - β (y), \mbox f or a r bi t r a r y (x, y) \in R^{u (N)} \times R^{s (N)},

Q (b_{c} (x) - b_{c} (y)) > 0 \mbox a s l o n g a s x \neq = y .

Q (b_{c} (x) - b_{c} (y)) > 0 \mbox a s l o n g a s x \neq = y .

W_{Z}^{s} (z_{0}, φ)

W_{Z}^{s} (z_{0}, φ)

W_{Z}^{u} (z_{0}, φ)

\frac{d u}{d t} = F (u),

\frac{d u}{d t} = F (u),

\frac{d a _{k}}{d t} = F_{k} (a) = (- \frac{k ^{4} π ^{4}}{L ^{4}} + \frac{λ k ^{2} π ^{2}}{L ^{2}} - λσ) a_{k} - \frac{λ k ^{2} π ^{2}}{L ^{2}} j_{1}, j_{2}, j_{3} \in Z j_{1} + j_{2} + j_{3} = k \sum a_{j_{1}} a_{j_{2}} a_{j_{3}} for all k > 0,

\frac{d a _{k}}{d t} = F_{k} (a) = (- \frac{k ^{4} π ^{4}}{L ^{4}} + \frac{λ k ^{2} π ^{2}}{L ^{2}} - λσ) a_{k} - \frac{λ k ^{2} π ^{2}}{L ^{2}} j_{1}, j_{2}, j_{3} \in Z j_{1} + j_{2} + j_{3} = k \sum a_{j_{1}} a_{j_{2}} a_{j_{3}} for all k > 0,

W \oplus T = {a \in H^{'} ∣ (a_{1}, \dots, a_{m}) \in W \mbox an d (a_{m + 1}, a_{m + 2}, \dots) \in T} \subset H^{'},

W \oplus T = {a \in H^{'} ∣ (a_{1}, \dots, a_{m}) \in W \mbox an d (a_{m + 1}, a_{m + 2}, \dots) \in T} \subset H^{'},

R^{\infty} \supset T = k > m \prod [a_{k}^{-}, a_{k}^{+}] .

R^{\infty} \supset T = k > m \prod [a_{k}^{-}, a_{k}^{+}] .

F_{k} (x) > 0 if x_{k} = a_{k}^{-}, \mbox an d F_{k} (x) < 0 if x_{k} = a_{k}^{+}, for all k > m .

F_{k} (x) > 0 if x_{k} = a_{k}^{-}, \mbox an d F_{k} (x) < 0 if x_{k} = a_{k}^{+}, for all k > m .

H^{'} = H_{1}^{'} \oplus H_{2}^{'} \oplus \dots,

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Code & Models

Repositories

https://bitbucket.org/dzako/dbcp_proof
noneOfficial

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

**Computer-assisted proof of heteroclinic connections

in the one-dimensional Ohta-Kawasaki model**

Jacek Cyranka*∗,‡,†, Thomas Wanner§*

‡ Institute of Computer Science and Computational Mathematics, Jagiellonian University

ul. S. Łojasiewicza 6, 30-348 Kraków, Poland

† Institute of Applied Mathematics and Mechanics, University of Warsaw

Banacha 2, 02-097 Warszawa, Poland

∗ Department of Mathematics, Rutgers, The State University of New Jersey, 110 Frelinghusen Rd, Piscataway, NJ 08854-8019, USA

§ Department of Mathematical Sciences, George Mason University, Fairfax, VA 22030, USA

[email protected], [email protected], [email protected]

Abstract. We present a computer-assisted proof of heteroclinic connections in the one-dimensional Ohta-Kawasaki model of diblock copolymers. The model is a fourth-order parabolic partial differential equation subject to homogeneous Neumann boundary conditions, which contains as a special case the celebrated Cahn-Hilliard equation. While the attractor structure of the latter model is completely understood for one-dimensional domains, the diblock copolymer extension exhibits considerably richer long-term dynamical behavior, which includes a high level of multistability. In this paper, we establish the existence of certain heteroclinic connections between the homogeneous equilibrium state and local and global energy minimizers.

The proof of the above statement is conceptually simple, and combines several techniques from some of the authors’ and Zgliczyński’s works. Central for the verification is the rigorous propagation of a piece of the unstable manifold of the homogeneous state with respect to time. This propagation has to lead to small interval bounds, while at the same time entering the basin of attraction of the stable fixed point. For interesting parameter values the global attractor exhibits a complicated equilibrium structure, and the dynamical equation is rather stiff. This leads to a time-consuming numerical propagation of error bounds, with many integration steps. This problem is addressed using an efficient algorithm for the rigorous integration of partial differential equations forward in time. The method is able to handle large integration times within a reasonable computational time frame, and this makes it possible to establish heteroclinic at various nontrivial parameter values.

Keywords:

Diblock copolymer model, dissipative partial differential equation, attractor structure, heteroclinic connection, multistability, computer-assisted proof.

AMS classification:

Primary: 35B40, 35B41, 35K55, 65C20. Secondary: 35Q99, 15B99, 68U20.

1 Introduction

The goal of this paper is to propose a computer assisted method of constructively proving existence of heteroclinic connections between stationary solutions of dissipative partial differential equations (PDEs). As a case study we establish a computer assisted proof of heteroclinic connections in the one-dimensional Ohta-Kawasaki diblock copolymer model. Precisely, we establish the existence of certain heteroclinic connections between the homogeneous equilibrium state, which represents a perfect copolymer mixture, and all local and global energy minimizers. In this way, we show that not every solution originating near the homogeneous state will converge to the global energy minimizer, but rather is trapped by a stable state with higher energy. This phenomenon can not be observed in the one-dimensional Cahn-Hillard equation, where generic solutions are attracted by a global minimizer.

Our method is general. To achieve the presented goal we develop a framework combining efficient algorithms and implementation in the C++ programming language performed by the first author [11] and the theory of cone conditions and self-consistent bounds developed by Zgliczyński et al. [2, 14, 15, 59, 60, 61, 63]. We remark that the software is not hard coded for the particular equation studied in this paper, but is rather a generic tool that can be easily adapted for study of any dissipative PDE on a one dimensional domain with a polynomial nonlinearity and periodic/Dirichlet/Neumann boundary conditions.

To the best of our knowledge, this is the first result about proving constructively the existence of heteroclinic connections for dissipative PDEs. In addition, it also rigorously shows for the first time that there are infinitely many orbits which are trapped by each of the two local and of the two suspected global energy minimizers. In fact, the proof of our result below shows that not all solutions of one dimensional Ohta-Kawasaki diblock copolymer model, which originate close to the homogeneous state, and satisfty certain space-translational symmetry, converge to the global energy minimizer.

Our approach differs in spirit from techniques based on Newton-Kantorovich type arguments in a functional analytic setting, which were used for example in [7, 20, 49]. We would like to mention that establishing the existence of connecting orbits in finite-dimensional continuous and discrete dynamical systems has been the subject of a number of studies, see for example [23, 49, 56, 57, 58]. In addition, the global dynamics of the Swift-Hohenberg model has been uncovered using a Morse decomposition based approach in [17], and the existence of periodic orbits for the ill-posed Boussinesq equation was established in [7, 15]. Finally, the existence of a heteroclinic orbit between the trivial and a nontrivial stationary state has been established for the one-dimensional $p$ -Laplace equation in [66]. Also, private communication revealed that Zgliczyński could verify the existence of a heteroclinic connection in the Kuramoto-Sivashinsky equation using similar techniques in the recent work in progress [64].

In the following two subsections, we first briefly introduce the Ohta-Kawasaki diblock copolymer model, structure of its dynamics, relevant bifurcation diagrams, and present our main result (Theorem 1.1). After that, we present the computer assisted framework used to prove the main result.

1.1 The diblock copolymer model description and the main result

Phase separation phenomena in materials sciences provide a rich source for pattern formation mechanisms. In many situations, the resulting models are dissipative parabolic partial differential equations which are amenable to rigorous mathematical treatment, and can therefore explain the underlying pattern formation principles. In the current paper, we consider one of these models which has been of significant interest in recent years. This equation is of phase field type, and it models microphase separation in diblock copolymers. From a mathematical perspective, the model has been described in [37], in its original form it was proposed by Ohta and Kawasaki [38], as well as Bahiana and Oono [1]. See also the discussion in [9, 10]. Consider a material constrained to the bounded domain $\Omega\subset\mathbb{R}^{d}$ . Then the model is based on an energy functional which is comprised of the standard van der Waals free energy with an additional nonlocal term, given by

[TABLE]

With this energy one can associate gradient-like dynamics, which is then used to model the evolution of the material as a function of time. The diblock copolymer model in its standard form considers the evolution equation

[TABLE]

where the nonlinearity is defined as $f(u)=-F^{\prime}(u)$ , and which in our situation is given by $f(u)=u-u^{3}$ . This evolution equation is based on the gradient in the $H^{-1}$ -topology and uses homogeneous Neumann boundary conditions for both $u$ and $\Delta u$ . The diblock copolymer model (1.1) is an extension of the celebrated Cahn-Hilliard equation [6], which corresponds to the special case $\sigma=0$ and serves as a fundamental model for the phase separation phenomena spinodal decomposition [31, 32, 44, 45, 52] and nucleation [3, 4, 18]. The quantity $\mu$ average of the solution over the domain is conserved in time.

For the present paper, we focus exclusively on the diblock copolymer model on one-dimensional domains $\Omega$ . Recent studies in [24] have demonstrated that the structure of its equilibrium set is considerably richer than the one for the Cahn-Hilliard equation described in [21]. This is visualized in Figure 1, which contains equilibrium bifurcation diagrams for the parameter values $\sigma=0$ , $3$ , $6$ , and $9$ , from top left to bottom right, and for total mass $\mu=0$ and domain $\Omega=(0,1)$ . Each diagram indicates the structure of the equilibrium set as a function of the new parameter $\lambda=1/\epsilon^{2}$ .

The interesting regime is the limit $\epsilon\to 0$ , and by switching to $\lambda$ the diagrams become easier to visualize. The horizontal straight lines in all images represent the constant function $u\equiv\mu=0$ , which is clearly an equilibrium, and which we call the trivial solution. As $\lambda$ increases, nontrivial solutions bifurcate from the trivial branch in pairs, even though we only plot one of these branches. This is due to the fact that the vertical axes of the bifurcation diagrams show the $L^{2}(\Omega)$ -norm of the stationary state, and the bifurcation pairs are related by a norm-preserving symmetry.

For the Cahn-Hilliard case $\sigma=0$ , i.e., in the top left bifurcation diagram, pairs of solutions are created at the bifurcation points $\lambda=k^{2}\pi^{2}$ for $k\in\mathbb{N}$ , and the solutions on the $k$ -th branch look qualitatively similar to the function $\cos(k\pi x)$ . In this case, each nontrivial solution branch corresponds to exactly two nontrivial solutions. It was shown in [21] that there are no other equilibrium solutions for the equation, i.e., none of the bifurcation branches exhibit secondary bifurcations. Moreover, only the solutions on the first branch are stable equilibria, and they correspond to solutions with one transition layer. The situation changes dramatically for $\sigma>0$ . In the diblock copolymer case, all branches which are created at the trivial solution exhibit secondary bifurcations. Even more importantly, these secondary bifurcations can lead to the creation of multiple stable equilibrium solutions. While as of yet there are no classical mathematical proofs for these statements, computer-assisted proof methods were developed in [55], which allow for the validation of these equilibrium branches. See also the work in [46, 54]. Some first results for two-dimensional domains can be found in the recent paper [51].

But how does this multistability manifest itself in the evolution equation (1.1)? Based on the underlying physical situation, one is usually interested in studying the long-term behavior of solutions of the diblock copolymer model which originate close to the (unstable) homogeneous equilibrium state $u\equiv\mu$ . These long-term limits are usually periodic in space with a specific period which depends on the precise parameter combinations. Since the system exhibits multiple stable states which all are at least local minimizers of the energy, it is natural to wonder whether typical solutions converge to the global energy minimizers or are trapped earlier. A systematic numerical study of this long-term behavior was performed in [24], and it was able to show that for fixed parameters $\lambda>0$ and $\sigma>0$ , most solutions starting near the homogeneous state usually converge to the same long-term limit. Moreover, the numerics indicated that this preferred long-term limit is usually not the stable equilibrium with the lowest energy. For more details, we refer the reader to [24]. In addition, the recent paper [53] contains a heuristic explanation of these numerical observations. We would like to point out, however, that to the best of our knowledge there is no rigorous complete description of the equilibrium set of the diblock copolymer model. Thus, even though the bifurcation diagrams in Figure 1 for $\sigma>0$ show multiple stable equilibria which can be ordered according to their energies, the lowest energy steady state is not known to be the global energy minimizer — there might be more equilibria than the ones shown in the diagram. However, since it is generally assumed that these diagrams are complete in the shown parameter ranges, we will refer in the following to the stable steady state in a diagram with the lowest energy as the suspected global minimizer.

Describing the complete attractor in a mathematically rigorous manner has proven to be elusive, even for the Cahn-Hilliard special case $\sigma=0$ . While the equilibrium structure of the Cahn-Hilliard equation is completely known for one-dimensional base domains [21], the precise structure of heteroclinic connections is an open question and has been answered only partially [4, 22]. Even for special two-dimensional domains such as the unit square, the structure of the Cahn-Hilliard attractor is extremely complicated, see for example the discussions in [19, 28, 29, 30]. In the case of the full diblock copolymer model, i.e., for $\sigma>0$ , even less is known. There are mathematical results on specific types of equilibrium solutions for two- and three-dimensional domains $\Omega$ , see for example [40, 41, 42, 43] as well as the references therein. Moreover, numerical studies have been performed of the long-term behavior of solutions of (1.1) which originate close to the homogeneous state [8, 9].

Despite the above numerical evidence, there has been no mathematical result which shows that large numbers of solutions of the diblock copolymer model which start near the trivial equilibrium solution are in fact trapped by a local energy minimizer, and that at the same time energy minimizers with lower energy can be reached from the homogeneous state as well. It is the goal of the present paper to close this gap for the bifurcation diagram shown in Figure 2.

This will be accomplished using a computer-assisted proof technique, and it will establish the existence of heteroclinic solutions between the homogeneous state and the local minimizers, as well as between the homogeneous state and local minimizers with lower energy, which are in fact suspected to be the global minimizers. More precisely, we verify heteroclinic connections from the state indicated by a yellow square, to both the equilibria indicated by yellow circles and the equilibria indicated by a yellow star. Our main result is the following theorem

Theorem 1.1 (Existence of heteroclinic connections).

Consider the diblock copolymer equation (1.1) on the one-dimensional domain $\Omega=(0,1)$ , for interaction lengths $\lambda=1/\epsilon^{2}=16\pi^{2}$ and $\sigma=16$ , and for total mass $\mu=0$ . Then there exist heteroclinic connections between the unstable homogeneous stationary state $u\equiv\mu$ and each of the two local energy minimizers which are indicated by yellow circles in Figure 2, and shown in the top right panel of Figure 3. In addition, there exist heteroclinic connections between the unstable homogeneous stationary state and each of the two suspected global energy minimizers which are indicated by yellow stars in Figure 2, and which are shown in the top left panel of Figure 3. In other words, for the above parameter values the diblock copolymer equation (1.1) exhibits multistability in the sense that local or suspected global energy minimizers can be reached from the homogeneous state.

We consider the specific parameter value $\lambda=16\pi^{2}\approx 157.91$ . It was shown in [24] that for this parameter value the diblock copolymer model (1.1) exhibits multistability which can be observed in practice. More precisely, if one randomly chooses initial conditions near the homogeneous state $u\equiv\mu=0$ , then some of the resulting solutions converge to the suspected global energy minimizer, while some are trapped in a local energy minimum At this value, the diblock copolymer model has eight nontrivial equilibrium solutions, which are shown in Figure 3.

The top left image shows the two suspected global energy minimizers, which are marked by a star in Figure 2, while the two local minimizers are indicated by a yellow circle, and are shown in the top right panel of Figure 3. Finally, the four colored solutions in the two lower panels of Figure 3 are of index one, and they all lie on the red secondary bifurcation branch in Figure 2. We would like to point out that the equilibrium solutions on the first and second branch bifurcating from the horizontal trivial solution line have two and three transition layers, respectively. The reason for this is described in more detail in [24]. The above theorem will be established as Theorem 4.1, which for computational reasons considers a rescaled version of the diblock copolymer model.

We remark that both of the local and global minimizers are symmetric. In Figure 3 it is clearly seen that the dominant harmonic of the local minimizer is $q=3$ , whereas the dominant harmonic of the global minimizer is $q=2$ (the $a_{q}\cos(qx)$ term of the solution’s expansion in Fourier series dominates). Also, the space of sequences $(a_{k})$ with $a_{k}=0$ for $k\neq q\mathbb{N}$ is invariant under the flow, which is also suggested by the numerical data from Appendix A. We make use of those symmetries in the construction of the stretched unstable isolating blocks, presented in Section 4. Precisely, to obtain the connection to the suspected global minimizer we stretch the block in the direction of the second eigenfunction $\cos{2x}$ , whereas to obtain the connection to the local minimizer we stretch the block in the direction of the third eigenfunction $\cos{3x}$ .

1.2 Framework for proving constructively connecting orbits in parabolic PDEs

In recent years, computer-assisted proof techniques have been used extensively in the context of nonlinear partial differential equations, and they are based on a number of different approaches. Here we present our approach for the computer assisted proof of Theorem 1.1.

Our method is general and can be used to achieve similar results for other dissipative PDEs. The crux of the whole algorithm for the proof of Theorem 1.1 is a procedure of rigorous propagation of a piece of the unstable manifold of the homogeneous state with respect to time. This propagation has to lead to small interval bounds, while at the same time entering the basin of attraction of the stable fixed point. For interesting parameter values the global attractor exhibits a complicated equilibrium structure, and the dynamical equation is rather stiff. This leads to a numerical propagation of error bounds involving lots of integration steps. For example the successful proof of the main result required performing of thousands of numerical integration steps. To address this problem we have developed an efficient algorithm for rigorous integration of dissipative PDEs forward in time [11]. Its highlight is an implementation of validated fast Fourier transforms, and generic and efficient implementation of a time-stepping scheme in C++ programming language combining several techniques, including the interval arithmetic, automatic differentiation, Taylor method, and Lohner’s algorithm [25, 26]. By calling it generic we mean that we developed a software tool that can be easily adapted to other dissipative PDEs, for example we have successfully integrated other equations including the Burgers, Swift-Hohenberg and Kuramoto-Sivashinsky PDEs [11, 12, 14]. Performing all required computations for the proof of the theorem presented in this paper takes only about 10 minutes on a laptop. All the numerical computations we performed are reproducible, the numerical data from the proofs, and the software codes with compile/use instructions are published online as a bitbucket repository [13]. We are convinced that the method is able to handle large integration times within a reasonable computational time frame, and in a future work we will explore applicability of our algorithm to parameter regimes in which the studied model exhibits much more complicated attractor structure and numerical stiffness. Explicit bounds for the unstable manifold and basins of attraction of the stationary solutions are computed using now standard techniques based on cone conditions [58, 62, 65] and logarithmic norms [12, 14, 59]. Computation of explicit bounds for the unstable manifold and basins of attraction of the stationary solutions is computationally very cheap compared to the numerical integration part.

1.3 Structure of the paper

The remainder of the paper is organized as follows. In Section 2 we describe the general setting for our approach, while Section 3 is devoted to its theoretical foundations. We recall basic definitions and results for isolating blocks, self-consistent bounds, as well as the cone condition. Based on this, Section 4 contains the proof of our main theorem. The next three sections are concerned with the computational aspects of our work. While the rigorous integration algorithm is described in Section 5, the following Section 6 shows how the cone condition can be verified in infinite dimensions. General remarks about the software can be found in Section 7. Section 8 contains conclusions and future plans. Finally, Appendix A.1 contains some numerical data from the computer assisted part of the proof. Throughout, from now on, we will drop the label “suspected” when we talk about the two suspected global minimizers in the above situation.

2 Problem formulation and basic setting

We begin by presenting the basic setup for our computer-assisted proof. As was mentioned in the introduction, our goal is to consider the diblock copolymer model (1.1) on the domain $\Omega=(0,1)$ . Unfortunately, for the parameter range of interest this choice of domain leads to extreme stiffness in the discretized equations. We will address this issue through rescaling techniques, and as a consequence our basic setup will be for a general one-dimensional domain of the form $\Omega=(0,L)$ . The length $L$ will be chosen later in such a way to reduce the stiffness of the problem. Thus, we study problem (1.1) in one dimension in the form

[TABLE]

which is equivalent to (1.1) up to a rescaling of time. In this formulation, we use the parameter $\lambda=1/\epsilon^{2}$ , and consider the partial differential equation subject to the mass and boundary constraints given by

[TABLE]

As before, the nonlinearity is given by $f(u)=u-u^{3}$ , and we consider the case $\mu=0$ of zero total mass, which means that we have equal amounts of the two polymers A and B.

Our basic discretization of this infinite-dimensional problem is based on the Fourier series expansion. More precisely, due to the imposed homogeneous Neumann boundary conditions we consider the expansion

[TABLE]

which is based on the cosine Fourier basis. We would like to point out that due to our choice $\mu=0$ for the total mass we automatically have

[TABLE]

Notice also that the domain length $L$ enters the expansion. We will see later that by choosing $L$ large enough the stiffness of the problem (3) can be reduced, since the eigenvalues of the linear part of a suitable finite-dimensional approximation are more averaged. To simplify notation, we will use the abbreviation

[TABLE]

for the basis functions in the above expansion. By default we assume that the coefficients $a_{k}$ are dependent on time $t$ , and we therefore frequently drop the explicit mentioning of the temporal argument. Upon substituting the series representation of $u$ into the cubic part of the nonlinearity $f$ one obtains

[TABLE]

The last expression can be further simplified using the fact that $\cos{\frac{\pi kx}{L}}=\frac{1}{2}\left(\exp\frac{i\pi kx}{L}+\exp\frac{-i\pi kx}{L}\right)\;$ , as long as we recall $a_{0}=0$ and define $a_{-k}=a_{k}$ for $k\in\mathbb{N}$ . By reordering one further obtains

[TABLE]

If we finally substitute this expression for $u^{3}$ , together with the original series expansion for $u$ , into the differential equation (3), then one obtains after division by $2$ and extraction of the coefficients of the basis functions $e_{k}(x)$ the identities

[TABLE]

In this way, we have reformulated the original parabolic partial differential equation as an infinite system of ordinary differential equations. This system will be used throughout the remainder of the paper. Recall one more time that the coefficient sequence $a_{k}$ satisfies

[TABLE]

To close this section, we present two central definitions. The first one is concerned with the function space which forms the basis of our computer-assisted proof.

Definition 2.1 (Sequence space with algebraic coefficient decay).

Let $H$ denote the space $\ell^{2}(\mathbb{Z},\mathbb{R})$ , i.e., elements $a\in H$ are sequences $a\colon\mathbb{Z}\to\mathbb{R}$ such that $\sum_{k\in\mathbb{Z}}{|a_{k}|^{2}}\leq\infty$ . In addition, let $\widetilde{H}\subset H$ denote the subspace of $H$ which is defined by

[TABLE]

where $\nmid x\nmid=|x|$ for $x\neq 0$ and $\nmid 0\nmid=1$ . Finally, let the space $H^{\prime}$ be given by

[TABLE]

Notice that all three spaces are Hilbert spaces.

From now on, all computations we perform in infinite dimensions will be constrained to the space $H^{\prime}$ . Needless to say, any actual numerical computations performed on a computer have to take place in finite dimensions, and for this we make use of Galerkin approximations of (5), which are the subject of the following final definition of this section.

Definition 2.2 (Galerkin approximation).

Consider the infinite system of ordinary differential equations given in (5), and let $n\in\mathbb{N}$ be arbitrary. Then the $n$ -th Galerkin approximation of (5) is defined as

[TABLE]

for any coefficient index $k=1,\dots,n$ . In other words, the $n$ -th Galerkin approximation is an $n$ -dimensional system of ordinary differential equations, which depends only on the coefficients $a_{1},\ldots,a_{n}$ .

3 Isolating blocks, self-consistent bounds, and the cone condition

In this section we briefly recall the theoretical framework that is used to establish the main result of this paper. More precisely, we review the concepts of cone conditions and self-consistent bounds as introduced by Zgliczyński. In Section 3.1, we present a finite-dimensional version of $h$ -sets, isolating blocks, and cone conditions. This is followed in Section 3.2 by a discussion of self-consistent bounds, which allows one to extend the notion of $h$ -sets, isolating blocks, and cone conditions to the Hilbert space setting. Finally, in Section 3.3 we present the concept of infinite-dimensional cone conditions, i.e. we show how to verify cone conditions on infinite-dimensional self-consistent bounds. Since the results presented in this section are known, our presentation will be short, albeit as self-contained as possible. For more detailed discussions of $h$ -sets, isolating blocks, and cone conditions we refer the reader to [58, 62, 65], and references therein.

The concepts introduced in the present section will be used in two different ways in the proof of our main theorem, which is the subject of the next Section 4. On the one hand, we use the fact that the construction of an unstable infinite-dimensional isolating block which satisfies a cone condition provides us with explicit bounds for the unstable manifold of the unstable steady state contained in the block. On the other hand, the construction of a stable infinite-dimensional isolating block which satisfies a cone condition guarantees the local uniqueness of the stable stationary solution contained in the block, as well as explicit bounds for the size of its basin of attraction.

3.1 Isolating blocks and cone conditions in finite dimensions

We begin by reviewing the concepts of isolating blocks and cone conditions in a finite-dimensional setting. For this, let $n\in\mathbb{N}$ be arbitrary. Then for any vector $z\in\mathbb{R}^{n}$ we let $\|z\|$ denote a norm on $\mathbb{R}^{n}$ , which does not have to be the standard Euclidean norm. The basic types of sets necessary for our discussion are interval sets. For this, consider real numbers $x_{k}^{\pm}$ for $k=1,\ldots,n$ . Then we define the interval set $[x]$ by

[TABLE]

For any point $x_{0}\in\mathbb{R}^{n}$ and arbitrary $r>0$ , we let $B_{n}(x_{0},r)=\{z\in\mathbb{R}^{n}\ |\ \|z-x_{0}\|<r\}$ denote the open ball of radius $r$ centered at $x_{0}$ , and we use the abbreviation $B_{n}=B_{n}(0,1)$ for the unit ball centered at the origin. For a linear map $A\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$ , we let ${Sp}(A)$ denote the spectrum of $A$ . When considering product spaces such as $\mathbb{R}^{u}\times\mathbb{R}^{s}$ , we write elements $z\in\mathbb{R}^{u}\times\mathbb{R}^{s}$ in the form $z=(x,y)$ , where $x\in\mathbb{R}^{u}$ and $y\in\mathbb{R}^{s}$ denote the first and second component vectors of $z$ , respectively, and the superscripts $s$ , and $u$ refer to the stable and unstable spaces respectively. Moreover, if $f\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$ denotes a continuously differentiable function, and if $Z\subset\mathbb{R}^{n}$ is a compact set, then we define

[TABLE]

Consider now an ordinary differential equation of the form

[TABLE]

Then we denote by $\varphi(t,z_{0})$ the solution of (7) which satisfies the initial condition $z(0)=z_{0}$ . Notice that since we assume the twice continuous differentiability of $f$ , the solution of this initial value problem is in fact uniquely determined.

After these preparations, we now turn our attention to the concept of $h$ -sets. Informally speaking, an $h$ -set is a subset of Euclidean space, which becomes a cube after a suitable choice of coordinate system. In addition, the letter $h$ alludes to the fact that this choice of coordinate system will lead to a phase portrait for the underlying ordinary differential equation which exhibits typical hyperbolic-like form. The concept of $h$ -sets originated in the work on covering relations for multi-dimensional dynamical systems, see [65]. From the perspective of our specific application, the considered $h$ -sets have a trivial structure in the sense that the cube is always expressed with respect to canonical coordinates. Consequently, also the cone conditions which will be introduced below are expressed in the canonical coordinate system. More precisely, in our application to the unstable homogeneous state, two of the dimensions will be distinguished as unstable, while the homeomorphism $c_{N}$ appearing in the following definition of an $h$ -set is just a scaling. Despite this fact, we now recall the appropriate definitions in their original form.

Definition 3.1 ( $h$ -set [65, Definition 1]).

An $h$ -set is a quadruple of the form $(N,u(N),s(N),c_{N})$ which satisfies the following three conditions:

•

The set $N$ is a compact subset of $\mathbb{R}^{n}$ ,

•

the natural numbers $u(N),s(N)\in\mathbb{N}_{0}$ satisfy the identity $u(N)+s(N)=n$ , and

•

the mapping $c_{N}\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$ , where we write $\mathbb{R}^{n}=\mathbb{R}^{u(N)}\times\mathbb{R}^{s(N)}$ , is a homeomorphism which satisfies

[TABLE]

In the following, we will slightly abuse notation and refer to an $h$ -set as just $N$ , rather than always listing the full quadruple. In addition, with each $h$ -set we associate the sets

[TABLE]

and we define $\dim(N):=n$ .

The above definition shows that an $h$ -set $N$ is the product of two closed unit balls in some appropriate coordinate system. The two integers $u(N)$ and $s(N)$ are called the nominally unstable and nominally stable dimensions, respectively. The subscript $c$ refers to the new coordinates given by the homeomorphism $c_{N}$ . Note that if we have $u(N)=0$ , then one immediately obtains $N^{-}=\emptyset$ , and the identity $s(N)=0$ implies $N^{+}=\emptyset$ . The next definition introduces the concept of a horizontal disk.

Definition 3.2 (Horizontal disk [58, Definition 10]).

Let $N$ denote an arbitrary $h$ -set. In addition, consider a continuous mapping $b\colon\overline{B_{u(N)}}\to N$ , and define $b_{c}=c_{N}\circ b$ . Then we say that $b$ is a horizontal disk in $N$ , if there exists a homotopy $h\colon[0,1]\times\overline{B_{u(N)}}\to N_{c}$ such that

[TABLE]

In order to handle the hyperbolic structure of $h$ -sets, we use the notion of cones, which are defined through a quadratic form defined on the $h$ -set in the following way.

Definition 3.3 ( $h$ -set with cones [65, Definition 8]).

Let $N\subset\mathbb{R}^{n}$ denote an arbitrary $h$ -set, and consider the associated splitting $\mathbb{R}^{n}=\mathbb{R}^{u(N)}\times\mathbb{R}^{s(N)}$ . Furthermore, let $Q\colon\mathbb{R}^{u(N)}\times\mathbb{R}^{s(N)}\to\mathbb{R}$ be a quadratic form which is given as the difference

[TABLE]

where $\alpha\colon\mathbb{R}^{u(N)}\to\mathbb{R}$ and $\beta\colon\mathbb{R}^{s(N)}\to\mathbb{R}$ are positive definite quadratic forms. Then the pair $(N,Q)$ is called an $h$ -set with cones.

Definition 3.4 (Cone condition [65, Definition 10]).

Consider an $h$ -set with cones given by $(N,Q)$ , as introduced in the previous definition. In addition, let $b\colon\overline{B_{u(N)}}\to N$ denote a horizontal disk as in Definition 3.2. Then we say that $b$ satisfies the cone condition (with respect to $Q$ ), if and only if for arbitrary points $x,y\in\overline{B_{u(N)}}$ we have

[TABLE]

Recall that we have $b_{c}\colon\mathbb{R}^{n}\to\mathbb{R}^{u(N)}\times\mathbb{R}^{s(N)}$ , and $b_{c}=c_{N}\circ b$ .

The remaining definitions of this section are concerned with standard dynamical systems notions. We begin by recalling the concept of isolating block from Conley index theory [34].

Definition 3.5 (Isolating block).

Let $N$ denote an arbitrary $h$ -set. We say that $N$ is an isolating block for the vector field $f:\mathbb{R}^{n}\to\mathbb{R}^{n}$ , if and only if $c_{N}$ is a diffeomorphism, if the sets $N^{+}$ and $N^{-}$ are local sections for $f$ , i.e., the vector field $f$ is transversal to $N^{\pm}$ , and if the following two conditions are satisfied:

•

For all points $x\in N^{-}$ there exists a constant $\delta>0$ such that $\varphi(t,x)\not\in N$ for all $t\in(0,\delta]$ .

•

For all points $x\in N^{+}$ there exists a constant $\delta>0$ such that $\varphi(t,x)\not\in N$ for all $t\in[-\delta,0)$ .

In other words, the closed set $N^{-}$ consists of all exit points for the flow associated with (7), and the closed set $N^{+}$ is the set of entry points.

Lemma 3.6 (Graph representation [62, Lemma 5]).

Let $(N,Q)$ denote an arbitrary $h$ -set with cones, and let $b\colon\overline{B_{u(N)}}\to N$ be a horizontal disk which satisfies the cone condition as in Definition 3.4. Then there exists a Lipschitz continuous function $y\colon\overline{B_{u(N)}}\to\overline{B_{s(N)}}$ such that $b_{c}(x)=(x,y(x))$ .

Definition 3.7 (Hyperbolic fixed point).

Let $z_{0}\in\mathbb{R}^{n}$ and consider again the ordinary differential equation (7). Then $z_{0}$ is called a hyperbolic fixed point for (7) if $f(z_{0})=0$ and if ${Re}{\lambda}\neq 0$ for all eigenvalues $\lambda\in{Sp}(df(z_{0}))$ , where $df(z_{0})$ is the Jacobian matrix of the vector field $f$ at the point $z_{0}$ , and ${Re}{\lambda}$ denotes the real part of the eigenvalue $\lambda$ .

Now consider a hyperbolic fixed point $z_{0}$ which is contained in a given set $Z\subset\mathbb{R}^{n}$ . Then we define the stable and unstable manifolds of $z_{0}$ as

[TABLE]

respectively.

Important for our later applications will be the following theorem. It shows that for hyperbolic fixed points, one can always find an $h$ -set $N$ with cones in such a way that the unstable manifold $W^{u}_{N}(z_{0})$ is a horizontal disk in $N$ , and that the stable manifold $W^{s}_{N}(z_{0})$ is a vertical disk.

Theorem 3.8 (Invariant manifolds as disks [62, Theorem 26]).

Suppose that the point $z_{0}$ is a hyperbolic fixed point for the ordinary differential equation (7). Moreover, let $Z\subset\mathbb{R}^{n}$ denote an open set with $z_{0}\in Z$ . Then there exists an $h$ -set $N$ with cones such that the following hold:

•

We have both $z_{0}\in N$ and $N\subset Z$ .

•

The $h$ -set $N$ is an isolating block for the vector field $f$ .

•

The unstable manifold $W^{u}_{N}(z_{0})$ is a horizontal disk in $N$ which satisfies the cone condition.

•

The stable manifold $W^{s}_{N}(z_{0})$ is a vertical disk in $N$ which also satisfies the cone condition.

Finally, the unstable manifold $W^{u}_{N}(z_{0})$ can be represented as the graph of a Lipschitz continuous function over the unstable space of the linearization of $f$ at $z_{0}$ , and this graph is tangent to the unstable space at $z_{0}$ . An analogous statement holds for the stable manifold $W_{N}^{s}(z_{0})$ .

In our application below we are interested in a slightly different situation. In our case, the existence of the fixed point is not known ahead of time. Its existence is, however, a consequence of the existence of an isolating block. Such a theorem in the context of maps has been proven in [62, Theorem 10]. In the context of differential equations it has been first established in [33], see also [28, 61].

3.2 Self-consistent bounds

We now turn our attention to the concept of self-consistent bounds, which have been used extensively in recent years [14, 15, 28, 59, 60, 61, 63]. In the following, we recall the basic definitions and fundamental ideas behind this concept, and describe an application which will be used for our main result.

In its most general form, the method of self-consistent bounds applies to arbitrary dissipative evolution equations in a suitable subspace of the sequence space $\ell^{2}$ , subject to certain admissibility conditions which depend on the specific notion of dissipativity. In other words, one can suppose that the underlying evolution equation takes the abstract form

[TABLE]

where $F$ is defined on a suitable subspace of $\ell^{2}$ and takes values in the sequence space. Rather than presenting the concept of self-consistent bounds in this general setting, we right away specialize to the situation of the present paper. Our base space is the Hilbert space $H^{\prime}\subset\ell^{2}$ which was introduced in Definition 2.1, and we consider the infinite system of ordinary differential equations given by (5). In other words, in our situation the abstract equation (8) can be written in coordinate form as

[TABLE]

where $a_{0}(t)\equiv 0$ . Without going into details, we note that the right-hand side is well-defined for all sequences $a\in H^{\prime}$ . In fact, the algebraic decay imposed through the subspace $\widetilde{H}$ in Definition 2.1 guarantees the four-time continuous differentiability of the function $u$ given by (4), and therefore standard regularity results for parabolic partial differential equations ensure the well-posedness of (9) in the space $H^{\prime}$ .

The basic idea behind the method of self-consistent bounds is the construction of a suitable subset of the phase space of (8), in our case the Hilbert space $H^{\prime}\subset\ell^{2}$ , for which one can easily determine the flow across its boundary. Usually, such a set is constructed as an infinite product of intervals for each of the coefficients $a_{k}$ of the sequences $a\in H^{\prime}$ , where, in order to be able to handle this set computationally, a bound for infinite number of them is provided by an algebraic decay. More precisely, let $a_{k}^{-}$ and $a_{k}^{+}$ denote the lower and upper bounds, respectively, for the $k$ -th coefficient $a_{k}$ , i.e., we will assume that elements of the constructed set satisfy $a_{k}\in[a_{k}^{-},a_{k}^{+}]$ , as long as $k$ is sufficiently large. In detail, self-consistent bounds are defined as follows.

Definition 3.9 (Self-consistent bounds).

Consider a set of the form

[TABLE]

which is the product of a finite-dimensional set $W\subset\mathbb{R}^{m}$ and an infinite-dimensional tail $T\subset\mathbb{R}^{\infty}$ . Then the set $W\oplus T$ forms self-consistent bounds for the evolution equation (9) if and only if the following conditions are satisfied.

•

The infinite-dimensional tail $T$ is given by an infinite product of intervals in the form

[TABLE]

Moreover, there exists a usually large integer $M>m$ such that $a_{k}^{-}\leq 0\leq a_{k}^{+}$ for all $k\geq M$ , i.e., all but finitely many of the intervals in $T$ contain zero.

•

There exists a constant $C>0$ and an exponent $s\geq 6$ such that $|a_{k}^{\pm}|\leq C/|k|^{s}$ for all $k>m$ .

•

The vector field $F$ points inwards on the parts of the boundary of $W\oplus T$ which are associated with the components in the tail $T$ in the sense that for all sequences $x\in W\oplus T$ we have

[TABLE]

The evolution equation (9) represents a dissipative partial differential equation in the sense discussed in [59, 60, 63]. Moreover, one can show that the right-hand side operator in (9) is in fact continuous on the set $W\oplus T$ . The block decomposition of the underlying Hilbert space $H^{\prime}$ is given by the infinite sum

[TABLE]

where the subspace $H^{\prime}_{k}$ is the orthogonal projection onto the subspace of $H^{\prime}$ which corresponds to the $k$ -th cosine basis vector via (4). For all of our considerations we will not use the standard induced norm on $H^{\prime}$ , but rather the so-called block-infinity norm given by

[TABLE]

This norm was introduced in [59]. For later applications, we also define the projection $P_{n}\colon H^{\prime}\to\mathbb{R}^{n}$ which maps $a$ to $(a_{1},\ldots,a_{n})$ for all $n\in\mathbb{N}$ . Finally, within the set of self-consistent bounds $W\oplus T\subset H^{\prime}$ for (9) we have

•

uniform convergence and existence of a solution for the infinite-dimensional system (9) in the sense that the orbits that stay in the set during a finite time interval converge in the block-infitnity norm,

•

the solutions of (9) are uniquely determined, and

•

one can derive a bound for the Lipschitz constant of the associated flow.

For more details on the justifications for each of these statements we refer the reader to [59].

3.3 Cone conditions in infinite dimensions

In this section we briefly describe how the finite-dimensional cone conditions from Section 3.1 can be interpreted in the infinite-dimensional self-consistent bounds setting. The involved results have appeared in two recent works on computer-assisted proofs for partial differential equations. On the one hand, they were used to establish a heteroclinic connection for the Kuramoto-Sivashinsky partial differential equation in [64]. In addition, in the recent thesis [36] they are fundamental for the computation of rigorous enclosures of unstable manifolds in the Cahn-Hillard equation.

As a first step, we have to extend the notion of $h$ -set to the infinite-dimensional self-consistent bounds setting. This can be accomplished as follows.

Definition 3.10 (Self-consistent bounds and $h$ -sets).

Suppose that $W\oplus T\subset H^{\prime}$ establishes self-consistent bounds for the evolution equation (9). Then the set $W\oplus T$ is an $h$ -set, if and only if the finite-dimensional part $W$ is an $h$ -set in the sense of Definition 3.1.

Our primary interest is the verification of cone conditions for self-consistent bounds for the evolution equation (9) for the purposes of establishing the existence of an equilibrium solution in the self-consistent bound set, and to recognize the unstable manifold at this stationary state as a suitable horizontal disk. This is the subject of the next theorem, which states that if the cone condition has been verified for self-consistent bounds $W\oplus T$ and a specific matrix $Q$ , then there exists a stationary point $z_{0}\in W\oplus T$ for (9), and the unstable manifold $W^{u}_{N}(z_{0})$ is a horizontal disk in $W\oplus T$ which satisfies the cone condition.

Theorem 3.11 (Consequences of the cone condition).

Consider self-consistent bounds $W\oplus T\subset H^{\prime}$ for (9), and assume that the orthogonal projection $P_{n}(W\oplus T)$ is an isolating block for the $n$ -th Galerkin approximation (6) of (9) for all $n>M$ . In particular, there is a finite number of unstable directions, all of them are confined to the finite dimensional part of $W\oplus T$ , and all the directions corresponding to the infinite dimensional tail are stable. In addition, let $Q$ denote an infinite diagonal matrix with $Q_{jj}=1$ for all indices $j$ corresponding to the unstable directions and $Q_{jj}=-1$ for the stable ones. Using this matrix, we assume that

[TABLE]

for all $w\neq 0$ of the form $w=x-y$ , where $x,y\in W$ .

Then there exists an equilibrium $z_{0}\in W\oplus T$ for (9) and the unstable manifold $W^{u}_{N}(z_{0})$ is a horizontal disk in $W\oplus T$ which satisfies the cone condition, where we set $N=W\oplus T$ . We use here an extended notion of $h$ -set as a product of a compact subset of $\mathbb{R}^{M}$ and an infinite tail. Moreover, the manifold $W^{u}_{N}(z_{0})$ can be represented as the graph of a Lipschitz continuous function over the unstable space.

For a detailed proof we refer the reader to [36, 64]. Basically, this theorem follows from the finite-dimensional version given in Theorem 3.8, combined with the fact that within self-consistent bounds $W\oplus T\subset H^{\prime}$ both the Fourier series of all terms appearing in the right-hand side $F$ of (9) and all partial derivatives in $dF$ converge uniformly. This allows one to establish the existence of unstable manifolds for all Galerkin approximations of (9) with uniform bounds, which in addition are contained in a compact set. One can now employ a standard convergent subsequence argument to establish the existence of an unstable manifold for the infinite-dimensional system (9). Moreover, the verification of positive definiteness of the infinite-dimensional matrix $[d\,P_{m}F(W)]^{T}Q+Q[d\,P_{m}F(W)]$ reduces to the verification of positivity of a finite number of inequalities, due to the polynomial growth of the diagonal which stems from the quartic increase of the eigenvalues of the linear part of (9). Finally, since the equation is dissipative, for sufficiently large $k$ the eigenvalues of the linear part of (9) are all negative, and therefore we have $Q_{kk}=-1$ for all these indices. In order to verify (11), we make use of the following lemma.

Lemma 3.12 (Verifying the cone condition).

If for some $\varepsilon>0$ and for all $k\in\mathbb{N}$ one can show that

[TABLE]

is satisfied, then the inequality (11) holds.

It will be the subject of Section 6 below how this infinite-dimensional condition can be verified computationally, and this will provide the unstable equilibrium which forms the initial point for our heteroclinic solutions.

Finally, we would like to mention that the existence of a stable steady state of (9) can be accomplished by checking suitable estimates for logarithmic norms, which were introduced in [16, 27]. In addition, such logarithmic norm estimates will provide explicit bounds for the size of the basin of attraction of the stable steady state. For this, we basically apply the theory, algorithms, and code from [12, 14, 59], and the interested reader is referred to these papers for more details. We do not further describe the details this verification procedure, since the logarithmic norm condition is equivalent to the cone condition (12) with only stable directions, i.e., with the matrix $Q$ being given by $-\text{Id}$ , where Id denotes the infinite-dimensional identity matrix.

4 Proof of the main result

With this section we begin proving our main result, which was described in Theorem 1.1 in the introduction. While this result guarantees the existence of heteroclinic connections between an unstable index two stationary solution, which in fact is the homogeneous steady state, and two stable local energy minimizers, its formulation concentrated on the diblock copolymer model on the domain $\Omega=(0,1)$ . This choice was motivated by the studies in [24], yet for numerical reasons it is not optimal. In fact, by naively choosing the parameters as in Theorem 1.1 one ends up with an extremely stiff Galerkin approximation (6).

In order to overcome this stiffness problem, we make use of standard rescaling arguments. For a fixed approximation dimension taking larger domain size results in ’averaging’ the set of linear part eigenvalues appearing in the approximation of a fixed dimension. For this, assume that $u$ denotes a solution of the diblock copolymer model (1.1) on the one-dimensional domain $\Omega=(0,1)$ , and for fixed parameters $\lambda=1/\epsilon^{2}>0$ and $\sigma>0$ . Furthermore, let $L>0$ be a fixed interval length, and define a function $v:\mathbb{R}_{0}^{+}\times(0,L)\to\mathbb{R}$ via

[TABLE]

Then one can easily verify that $v$ solves the diblock copolymer equation

[TABLE]

which is exactly of the form given in (3), but with rescaled parameters $\lambda$ and $\sigma$ .

It turns out that for the parameter combinations given in Theorem 1.1 it is most convenient to choose the new domain length $L=2\pi$ , since this reduces the size of the eigenvalues of the linearized Galerkin approximation (6). Thus, in the following, we will establish the following rescaled result.

Theorem 4.1 (Rescaled version of the main theorem).

Consider the diblock copolymer equation (3) on the domain $\Omega=(0,L)$ with $L=2\pi$ , for interaction lengths $\lambda=4$ and $\sigma=4/\pi^{2}\approx 0.405284735$ , and for total mass $\mu=0$ . Then there exist heteroclinic connections between the unstable homogeneous stationary state $u\equiv 0$ , which is indicated by a yellow square in Figure 4, and

(a)

Each of the two local energy minimizers which are indicated by yellow circles, as well as 2. (b)

Each of the two global energy minimizers which are indicated by yellow stars.

In other words, for the above parameter combinations, the diblock copolymer equation exhibits multistability, and all energy minimizers can be reached from any arbitrarily small neighborhood of the homogeneous state.

To illustrate the formulation of the above theorem, we have included the bifurcation diagrams for the new, rescaled situation in Figure 4. Compare also their original counterparts in Figure 2.

Before turning our attention to the proof of Theorem 4.1, we need to introduce some notation. In the following, we denote the unstable homogeneous steady state of (3) at origin by $z_{0}^{u}\in H^{\prime}$ . This equilibrium is indicated by a yellow square in Figure 4. In addition, let $z_{0}^{s,loc}\in H^{\prime}$ denote one of the two stable local energy minimizers of (1), which are indicated by a yellow circle in Figure 4. Let $z_{0}^{s,glob}\in H^{\prime}$ denote one of the two global energy minimizers of (1), which are indicated by yellow star in Figure 4. One can easily see that the theorem only has to be verified for one of these local/global minima, the result for the second one follows from symmetry arguments.

For each of the involved equilibrium solutions, we need to construct self-consistent bounds containing them. Thus, we let

[TABLE]

denote an infinite-dimensional stable isolating block which contains the stable steady state $z_{0}^{s,loc}$ , and which forms self-consistent bounds. We also let

[TABLE]

denote the infinite-dimensional stable isolating block which contains the stable steady state $z_{0}^{s,glob}$

As long as it also satisfies the cone condition, the existence of stable steady states in $\mathcal{W}_{\text{s}}^{loc}$ and $\mathcal{W}_{\text{s}}^{glob}$ follows then from Theorem 3.11 with all directions being stable directions, as illustrated in Figure 5. Now consider the unstable equilibrium involved in the heteroclinic. Let

[TABLE]

denote the infinite-dimensional unstable isolating block which contains the known unstable equilibrium $z_{0}^{u}=0$ , and which forms self-consistent bounds. In this case, if the block satisfies the cone condition, then the existence of a two-dimensional unstable manifold is guaranteed at the origin by Theorem 3.11, as is the fact that the manifold is a horizontal disk in $\mathcal{W}_{\text{u}}$ . See also the illustration in Figure 6. After these preparations, we can now proceed with the proof of Theorem 4.1.

Proof.

The verification of Theorem 4.1 is divided into three parts, which are described in more detail below. We would like to point out already now that all sets used in the proof are constructed using rigorous numerical programs which are based on interval arithmetic, see for example [35]. In this way, all computations done by a computer are validated in a strict mathematical sense. In particular, if the program verifies the cone condition for the whole set under consideration, this automatically implies that the condition is satisfied for all elements of this set. For more illustrations of this approach, also in other contexts, we refer the interested reader to [48].

First step: Using a computer program based on interval arithmetic, we first construct self-consistent bounds $\mathcal{W}_{\text{u}}\subset H^{\prime}$ in such a way that $\mathcal{W}_{\text{u}}$ is an unstable isolating block for (9). For this, the set $\mathcal{W}_{\text{u}}$ is centered at the origin, which represents our known unstable equilibrium $z_{0}^{u}$ for (9). More specific information regarding the specific choice of $\mathcal{W}_{\text{u}}$ is presented in Appendix A.1, which contains verified numerical data from the program. These results show that the topological structure of $\mathcal{W}_{\text{u}}$ is exactly of the form shown in Figure 6. In addition, the set $\mathcal{W}_{\text{u}}$ is verified to satisfy the cone condition (12), using the method presented in Section 6. Now Theorem 3.11 implies that the unstable manifold $W_{N}^{u}(z_{0}^{u})$ is a horizontal disk in $\mathcal{W}_{\text{u}}$ , and therefore the unstable manifold at the origin crosses one of the faces of $\mathcal{W}_{\text{u}}$ . Let $W_{0}$ denote such a face, see also the diagram in Figure 7. The whole construction of $\mathcal{W}_{\text{u}}$ relies on the fact that the unstable equilibrium $z_{0}^{u}$ is given by the origin. We can therefore determine the eigenfunctions of the linearization of (9) at $z_{0}^{u}$ explicitly, and this easily shows that all eigenfunctions are just the basis functions $e_{k}(x)$ for $k\in\mathbb{N}$ . Furthermore, recall from our brief discussion in the introduction that for our choice of $\sigma$ the two unstable eigendirections correspond to the eigenfunctions $e_{2}(x)$ and $e_{3}(x)$ . In order to make the subsequent parts of the proof computationally faster, the set $\mathcal{W}_{\text{u}}$ is chosen in such a way that it is much more elongated in the second eigenfunction direction compared to all the other directions, which is given by $e_{3}(x)$ (we denote the face of $\mathcal{W}_{\text{u}}$ in this direction by $W_{0}(e_{3})$ ) for the case of the connection to the local minimizer. Additionally, the set $\mathcal{W}_{\text{u}}$ is chosen in such a way that it is much more elongated in the third eigenfunction $e_{2}(x)$ direction for the case of the connection to the global minimizer (we denote the face of $\mathcal{W}_{\text{u}}$ in this direction by $W_{0}(e_{2})$ ). See again Figure 7.

Second step: Similar to the first step, we use interval arithmetic to construct self-consistent bounds $\mathcal{W}_{\text{s}}^{loc},\mathcal{W}_{\text{s}}^{glob}\subset H^{\prime}$ in such a way that $\mathcal{W}_{\text{s}}^{loc},\mathcal{W}_{\text{s}}^{glob}$ is a stable isolating block for (9). The block $\mathcal{W}_{\text{s}}^{loc}$ is centered at a numerically computed approximation to the actual stable equilibrium $z_{0}^{s,loc}$ , and the topological structure of $\mathcal{W}_{\text{s}}^{loc}$ is as indicated in Figure 5. Analogously, the block $\mathcal{W}_{\text{s}}^{glob}$ is centered at a numerically computed approximation to the stable equilibrium $z_{0}^{s,glob}$ – the global minimizer. The finite-dimensional part of this set is determined in an appropriate coordinate system, which is given by the numerically determined eigenbasis of the Jacobian matrix at the numerical approximation of $z_{0}^{s,loc}$ , $z_{0}^{s,glob}$ respectively. As a consequence, the block decomposition (10) of $H^{\prime}$ takes on the following form. The first $m$ coordinates correspond to the eigenbasis of the $m$ -th Galerkin approximation (6), while the remaining coordinate directions are given by the standard basis functions $e_{k}(x)$ for $k>m$ . In other words, we have to take into account that the Jacobian at $z_{0}^{s,loc}$ , $z_{0}^{s,glob}$ is considered in two different coordinate systems, and this is handled as outlined in more detail in [12, 59]. Next, the computer program is used to compute a verified upper bound on the logarithmic norm of the Jacobian in block coordinates (10). As can be seen from the numerical data in Appendix A.1, this bound turns out to be negative. Notice that determining an upper bound for a logarithmic norm is equivalent to computing the cone condition (12) in the situation of only stable directions, and our program makes use of the techniques described in more detail in [12, 59]. Thus, the numerical data presented in the appendix ensures the existence of an attracting equilibriums for (9), namely $z_{0}^{s,loc}$ in $\mathcal{W}_{\text{s}}^{loc}$ , and $z_{0}^{s,glob}$ in $\mathcal{W}_{\text{s}}^{glob}$ . In addition, we obtain that the complete sets $\mathcal{W}_{\text{s}}^{loc}$ , $\mathcal{W}_{\text{s}}^{glob}$ are contained in the basin of attraction of $z_{0}^{s,loc}$ , and $z_{0}^{s,glob}$ respectively.

Third step: Since the first two steps of the proof have established sets $W_{0}(e_{2}),W_{0}(e_{3})$ which contains part of the unstable manifold of $z_{0}^{u}$ , as well as a subsets $\mathcal{W}_{\text{s}}^{loc}$ , $\mathcal{W}_{\text{s}}^{glob}$ which are contained in the basin of attraction of $z_{0}^{s,loc}$ , $z_{0}^{s,glob}$ respectively, it only remains to show that all solutions of (9) which originate in $W_{0}(e_{3})$ end up in $\mathcal{W}_{\text{s}}^{loc}$ , and all solutions of (9) which originate in $W_{0}(e_{2})$ end up in $\mathcal{W}_{\text{s}}^{glob}$ . This is accomplished in the last step of the proof through rigorous forward integration of the flow associated with (9). For this, we use the algorithm for rigorous integration described in the series of works [12, 25, 60, 63]. In fact, we make use of an optimized version of this algorithm which was developed in [11], and this version will be described in more detail in the next Section 5. In this way, we are able to rigorously propagate the sets $W_{0}$ forward in time in such a way that all computational errors are accounted for, and we obtain that the flow of (9) in $H^{\prime}$ transports $W_{0}(e_{3})$ into $\mathcal{W}_{\text{s}}^{loc}$ , and $W_{0}(e_{2})$ into $\mathcal{W}_{\text{s}}^{glob}$ . A cartoon of this procedure is shown in Figure 8. In the figure, the sets $[\Phi(kh,W_{0})]$ denote rigorous enclosures for the images of $W_{0}$ under the flow at time $kh$ , for integers $k>0$ . After $N+1$ steps of this forward integration by time $h$ , the final enclosure is shown to be contained in $\mathcal{W}_{\text{s}}^{loc}$ , and therefore the basin of attraction of $z_{0}^{s,loc}$ . In the actual computation performed by a computer program any initial condition in the face $W_{0}(e_{3})$ of $\mathcal{W}_{\text{u}}$ at time $T=3.02$ (after performing $1510$ numerical integration steps) is verified to be within $\mathcal{W}_{\text{s}}^{loc}$ , and any initial condition in the face $W_{0}(e_{2})$ of $\mathcal{W}_{\text{u}}$ at time $T=4.62$ (after performing $2310$ numerical integration steps) is verified to be within $\mathcal{W}_{\text{s}}^{glob}$ . As a consequence, one obtains the existence of a manifold $W^{u}(z_{0}^{u})\cap W^{s}(z_{0}^{s,loc})$ which connects the unstable equilibrium $z_{0}^{u}$ with the stable equilibrium $z_{0}^{s,loc}$ . As well as the existence of a manifold $W^{u}(z_{0}^{u})\cap W^{s}(z_{0}^{s,glob})$ which connects the unstable equilibrium $z_{0}^{u}$ with the stable equilibrium $z_{0}^{s,glob}$ For more details on the numerical data from the actual computations we refer the reader to Appendix A.1 and the numerical data available in an online repository [13]. Needless to say, the third and last step of the proof is the most time consuming part. ∎

We would like to close this section with a brief comment on the applicability of the above approach. The proof given above defines an algorithm for proving Theorem 4.1. The algorithm is general. While in the proof we have considered a very specific instance of the diblock copolymer model, for a particular set of parameters, in principle the algorithm applies to the verification of any saddle to sink connection for the diblock copolymer model (9). The algorithm describes numerical verification needed to be performed with a clear separation of the computation into three consecutive steps. We have implemented a computer program, realizing this algorithm, which is part of the publication [13] and available online. The input to the program are the parameters in the equation, and if the program terminates successfully we have proved the existence of a heteroclinic connection. The main prerequisite for termination is of course that the set propagated by the rigorous numerical integrator is mapped completely inside the stable isolating block $\mathcal{W}_{\text{s}}^{loc}$ . However, the process of obtaining such a proof for arbitrary parameters is not fully automated. In general, one has to try to make the set $\mathcal{W}_{\text{s}}^{loc}$ as large as possible, while at the same time keeping $W_{0}$ small. In this way, the accumulated effect of roundoff errors and inherent expansion in the differential equation can often be controlled.

5 Rigorous integration algorithm

To rigorously integrate (9) forward in time we apply an algorithm based on the rigorous Taylor method, as well as the Lohner algorithm for an efficient representation of sets and to avoid the wrapping effect. The Lohner algorithm is presented in detail in [25, 26]. We specifically use the algorithm introduced in [11], which allows for the efficient implementation of the Taylor method for differential equations in combination with the variational equations arising in partial differential equations. This algorithm combines automatic differentiation with the fast Fourier transform algorithm. In the remainder of this section, we briefly explain the ideas behind this algorithm, and discuss why it is much more efficient than the direct method of computation of normalized derivatives. For a more detailed discussion, we refer the reader to [11], where a computational complexity analysis can be found in the context of partial differential equations with a quadratic nonlinearity. We will demonstrate below that this analysis can be extended to the diblock copolymer model (3) with a cubic nonlinearity. In fact, in the situation of the present paper, the method of [11] leads to even greater speedups. One of the main advantages of the method of [11] is that it allows one to use a high-order Taylor method as opposed to the direct integration method. In fact, in the actual proof of the results of this paper we use a Taylor method of order $16$ .

To describe our approach in more detail, consider again the right-hand side $F$ of (9). In addition, we denote the right-hand side of the $m$ -th Galerkin approximation (6) by $F_{m}$ . We begin by giving a brief description of how the actual solution of an initial value problem for (9) can be enclosed in an efficient way. This is achieved by bounding the solution set of the differential inclusion given by

[TABLE]

where $z(t)$ takes values in $\mathbb{R}^{m}$ . More precisely, the solution set of (13) is bounded by first enclosing the solution of the Cauchy problem associated with (13), and then adding appropriate a-posteriori bounds, see for example [12, 25] for more details.

In order to rigorously enclose the solution of the Cauchy problem associated with (13), we make use of the classical Taylor method. Thus, the enclosure for the solution is propagated by one time step through the knowledge of the normalized derivatives of orders up to $p$ , where the integer $p$ denotes the order of the Taylor method. These normalized derivatives are defined as follows.

Definition 5.1 (Normalized derivative).

Let $n>0$ be an arbitrary integer, and let $w\colon[0,t_{max})\to\mathbb{R}$ denote a sufficiently smooth function. Then the $n$ -th normalized derivative of $w$ is defined as

[TABLE]

where $w^{(n)}(t)$ denotes the classical $n$ -th derivative of $w$ .

According to this definition, the normalized derivatives are precisely the terms occurring in a standard Taylor expansion. They can be determined efficiently via a well-known recursive formula from automatic differentiation, which is recalled in the next remark.

Remark 5.2 (Recursion for normalized derivatives).

For an ordinary differential equation of the form

[TABLE]

the normalized derivatives of $z$ can be determined as follows. If we set $z^{[0]}=z^{(0)}$ , then for arbitrary differentiation orders $j=1,\dots,q-1$ the normalized derivative $z^{[j]}$ is given by the recursive formula

[TABLE]

as long as $G$ is $q$ -times differentiable.

After these preparations we can now analyze the cost of computing the $(p+1)$ -st order Taylor method for the finite system of ordinary differential equations given by the Galerkin approximation (6). If we assume that the normalized derivatives of $u$ of orders $1,2,\dots,p$ have already been computed, then the cost of determining the $(p+1)$ -st normalized derivative is equivalent to the cost of evaluating the convolutions

[TABLE]

where the function $u$ is given by its time-dependent Fourier coefficients, and $[p]$ denotes $p$ -th normalized derivative with respect to time.

Since we are interested in using the Lohner algorithm [26], we not only have to keep track of the Fourier coefficients $\{a_{k}\}_{k=1}^{m}$ of the solution $u$ , but also of the associated partial derivatives vector

[TABLE]

for each $j=1,\dots,m$ the first order multi-variate polynomial with respect to an auxiliary variable $\xi=\left(\xi_{1},\dots,\xi_{m}\right)$

[TABLE]

is usually referred to as the jet. Thus, in order to emphasize that we work with jets, not scalar values we will from now on write $\mathbf{u}$ instead of $u$ . An algebra of jets can be formed defining addition and multiplication as the usual operations on polynomials. But the terms of the second order in $\xi$ are truncated after the multiplication. Our goal is the computation of the convolutions

[TABLE]

which can be expanded as

[TABLE]

If we now let $N$ denote the number of coefficients that are calculated in (14), along with their partial derivatives for the variational equations, then one can easily see that we basically have

[TABLE]

up to a constant, which has to be large enough to avoid the aliasing effect (in our algorithm the constant is equal to $2$ ). With this, one can finally analyze the computational cost of evaluating (14). While we do not present the details, the result is collected in the following remark.

Remark 5.3 (Computational cost of evaluating the convolution).

The total cost of the direct evaluation of the convolution given in (14) is proportional to $p^{2}N^{3}$ , where $N=m^{2}+m$ up to a constant. An example analysis can be found in [11].

Besides the above direct computation, it is well-known that one can evaluate convolutions as in (14) indirectly through the use of the fast Fourier transform algorithm. If we assume again that all normalized derivatives of $u$ of order $1,2,\dots,p$ have already been computed, then the discrete functional values of the normalized derivatives of $u$ , which are also called the $L$ -coefficients in [11], have already been determined, and they can therefore be cached in memory. We denote these values by

[TABLE]

where $\mbox{FFT}_{j}\left(\cdot\right)$ denotes the $j$ -th component of the discrete Fourier transform computed by the fast Fourier transform method, and the constant $c$ is specified below. Using these values, one can then determine (14) by computing the convolution for the $L$ -coefficients in the form

[TABLE]

and then applying the inverse fast Fourier transform, which we denote by $\mbox{FFT}^{-1}$ , in order to obtain

[TABLE]

After this, we can also calculate the $L$ -coefficients of $\widehat{\mathbf{u}}^{[p+1]}$ and cache them in memory, provided we intend to increase the order further, i.e., we evaluate

[TABLE]

While at first glance the approach based on the fast Fourier transform seems to be more involved, its computational savings are significant. To see this, recall that the cost of performing a fast Fourier transform is of the order $N\log{N}$ , and therefore the same is true for the computation of (16). Moreover, the cost of computing the convolution for the $L$ -coefficients (15) is just the cost of multiplying $p^{2}cN$ values, and this leads to the following observation.

Remark 5.4 (Computational cost of the convolution via the fast Fourier transform).

If we assume that the $L$ -coefficients of $\widehat{\mathbf{u}}$ have already been computed and cached in memory, then the total cost of the fast Fourier transform approach to the computation of (14) is proportional to

[TABLE]

An example analysis can be found in [11].

Typically, the integer $N$ is significantly larger than the order $p$ , since it counts the number of equations and variationals. Clearly the cost provided in Remark 5.3, i.e. $p^{2}N^{3}$ dominates both of the terms appearing in (17) – the cost provided in Remark 5.4 for any $p$ and $N$ . Also, increasing the order $p$ in case of the cost given in Remark 5.4 does not result in as severe increase of the total computational cost, as in the case of Remark 5.3. The typical choice of the constant $c$ in the above discussion is $c=3$ . Using some optimizations that were described in [11], one can in fact reduce the constant to $c=2$ .

Needless to say, in practice the fast Fourier transform technique greatly reduces the total running time of our computer-assisted proof.

Remark 5.5.

The actual running time of the integration step of the proof with parameters as reported in Appendix A.1 is triple reduced. The real time of the integration procedure on a laptop is 3m55s FFT vs 11m55 direct (without FFT) in case of the connection to the local minimizer, and 6m34s FFT vs 17m29s in case of the connection to the global minimizer.

We remark that the actual proof was performed with a moderate Galerkin approximation dimension ( $m=15$ ), which allows to achieve the proof in a reasonable time on a laptop. In particular performing of an analogous proof, but for a parameter regime further from zero will require larger approximations, due to much more complicated dynamics on the attractor. The analysis of the algorithm presented in this section shows that for larger approximation dimension even a higher speedup is anticipated.

6 Verifying the cone condition in infinite dimensions

This section is devoted to illustrating how the cone condition (12) can be verified in infinite dimensions. For the sake of brevity, this will only be done in the context of the set $\mathcal{W}_{\text{u}}$ . As mentioned earlier, the verification of the cone condition for the remaining set $\mathcal{W}_{\text{s}}^{loc}$ basically amounts to obtaining logarithmic norm estimates, and this can be achieved directly with the algorithms from [12, 59], and is therefore omitted in the following. To begin with, we would like to remind the reader that according to (12), we have to verify the estimate

[TABLE]

for the unstable isolating block $\mathcal{W}_{\text{u}}$ , for some $\epsilon>0$ and for all $k\in\mathbb{N}$ .

In order to accomplish this, the block $\mathcal{W}_{\text{u}}\subset H^{\prime}$ is divided into two parts, and we refer to this representation again as self-consistent bounds. More precisely, we have that

•

the finite part for $k=1,\dots,M$ is represented as a finite-dimensional interval product, i.e., the projection $\left(\mathcal{W}_{\text{u}}\right)_{k}$ is an interval, and

•

the tail, which corresponds to $k=M+1,\dots$ , is represented by an algebraically decaying bound, i.e., for all $x\in\left(\mathcal{W}_{\text{u}}\right)_{k}$ the inequality $|x|\leq C/|k|^{6}$ is satisfied. In order to explicitly relate the constant $C$ to the block $\mathcal{W}_{\text{u}}$ , we will use the notation $C(\mathcal{W}_{\text{u}})$ in the following.

Notice that $\mathcal{W}_{\text{u}}$ forms self-consistent bounds with the fixed polynomial decay, therefore the tail in the algorithm is represented by only the value of the constant $C>0$ .

In order to establish the cone condition, we first determine self-consistent bounds $\mathcal{C}\subset H^{\prime}$ , which bound all possible convolutions of elements from the block $\mathcal{W}_{\text{u}}$ , i.e., such that for all $k=0,\dots,M$ we have

[TABLE]

Since the nonlinearity in (9) is cubic, one therefore immediately obtains an enclosure for the partial derivatives matrix $DF$ in the form

[TABLE]

Before proceeding, we need to introduce some notation. if $I\subset\mathbb{R}$ denotes an arbitrary interval, then we denote the maximal absolute value of all points from $I$ by $|I|_{\max}$ . We can now estimate the second expression in (18) by

[TABLE]

The right-hand side in this estimate can be bounded further. For this, we split it into two parts to obtain

[TABLE]

The first part of this sum consists of expressions which contain at least one term from the finite part. In addition, the first summand is a finite sum, and we can therefore compute it term by term using interval arithmetic. In contrast, the second summand is infinite and consists exclusively of terms from the tail part, which can therefore be estimated. Using standard results from rigorous numerics for partial differential equations, due to a high regularity of the Fourier coefficients sequences all of the appearing sums have a finite value, which can be estimated efficiently.

Finally this yields the following estimate

[TABLE]

whose right-hand side is given explicitly — and from this one can easily establish the validity of (18) for $k=1,\dots,M$ .

For $k>M$ using now standard arguments, see e.g. [63], it can be shown that (18) is decreasing with $k$ for $k>M$ , due to the fast decay of eigenvalues. Hence, the negativity of (18) for $k>M$ follows from the negativity of (18) for $k=M$ .

Finally the cone condition can be verified.

7 Software

The computer program performing all computational steps required to prove Theorem 4.1, as described in more detail in the respective proof in Section 4, has been written in the C++ programming language. The source codes, some compiled binary files, and the numerical output from the actual proofs are available in an online bitbucket repository [13]. The published source codes for this work utilizes some parts of the rigorous numerics software package developed by the first author, which was dedicated to the study of dissipative partial differential equations based on Fourier series expansions. This package has been used in a number of different contexts as well, see for example [11, 12, 14].

The software package is written in such a way that it should be easily adaptable to other situations as well. More precisely, by performing only small modifications in the code, it can be applied to a range of dissipative partial differential equations which are written in terms of infinite Fourier series. Examples include the viscous Burgers equation, the Kuramoto-Sivashinsky equation, and the Swift-Hohenberg equation. The main part of the software package is a dedicated algorithm for performing time integration using Taylor’s method, and the efficient computation of the (higher order) time derivatives of the flow. The latter are required by Taylor’s method, and they are determined using the FFT algorithm described earlier.

8 Conclusion and future work

In this paper, we have presented a first instance in which heterolinic solutions can be rigorously determined for the diblock copolymer model on one-dimensional domains. Moreover, our results have shown that the model exhibits strong multistability in the sense that all (local or global) energy minimizers can be reached from arbitrarily small neigborhoods of the homogeneous steady state. In the near future we plan to extend the presented research to prove heteroclinic connections for a wide range of other parameter regimes, with particular emphasis on regions further from the homogeneous steady state. Moreover, we plan to perform similar studies on simple two-dimensional domains. Another intriguing possibility for future work is to apply the parametrization method [5] in the first step of the proof of Theorem 4.1. Recently, the parametrization method has been generalized to allow for the validation of unstable manifolds for parabolic partial differential equations, see for example [39, 50]. Using the parametrization method, bounds for the unstable manifold are obtained by computing its expansion in the form of a finite Taylor polynomial. In contrast, in the present work, we obtain bounds for the unstable manifold by verifying the cone condition on a linear segment. For nonlinear partial differential equations and nontrivial parameter values, it is expected that the unstable manifold behaves linearly only very close to the unstable steady state. We are convinced that, using the parametrization method, especially for high Taylor orders, one can compute bounds for the unstable manifold of the diblock copolymer model that reach much further than those obtained using a linear block and cone conditions. This in turn is expected to greatly reduce the work needed for performing the final step of the proof, i.e., using the rigorous integrator to demonstrate that the unstable manifold enters the basin of attraction of the stable steady state. This should make it possible to then establish a result analogous to Theorem 4.1, but for large parameter values.

9 Acknowledgments

The presented work was performed while the first author held post-doctoral positions at the Warsaw Center of Mathematics and Computer Science, Poland, and subsequently at Rutgers, The State University of New Jersey, USA. The idea for the study arose during the workshop “Rigorous computation for infinite-dimensional nonlinear dynamics” at the American Institute of Mathematics in Palo Alto, California, in August 2013. J. Cyranka was supported in part by NSF grant DMS 1125174 and DARPA contract HR0011-16-2-0033. The work of the second author was partially supported by NSF grants DMS-1114923 and DMS-1407087.

Appendix A Numerical data from the computer proof

A.1 Heteroclinic to the local minimizer

In this appendix we present some numerical data from the computer-assisted part of the proof presented in this paper. For the sake of brevity we list only a small number of the relevant modes which are involved in the constructed infinite-dimensional bounds. In addition, only three decimals of all numerical values are shown. Complete and detailed numerical data presented with larger precision can be found online http://bitbucket.org/dzako/dbcp_proof.

Unstable isolating block $\mathcal{W}_{u}\subset H^{\prime}$

around zero.

This is the first step of the proof of the main theorem from Section 4. We mark with text in bold the numerical values for the eigendirection $e_{3}$ along which the unstable manifold of the homogeneous state connects to the stable local minimizer ( $\mathcal{W}_{u}$ is stretched along this direction)..

[TABLE]

The infinite-dimensional cone condition (18) is satisfied with

[TABLE]

The file unstable_box_log.txt which can be found online [13] contains more detailed numerical data related to this computation. The actual running time of the program computing this set is negligible compared to the last step, i.e., the numerical integration. Therefore, we do not report it.

Stable isolating block $\mathcal{W}_{\text{s}}^{loc}\subset H^{\prime}$

This is the second step of the proof of the main theorem from Section 4. The block $\mathcal{W}_{\text{s}}^{loc}$ is centered at a numerically computed approximation to the actual stable equilibrium $z_{0}^{s,loc}$ which was loaded from the provided file fixedPoint.in.

[TABLE]

The infinite-dimensional cone condition (18) is satisfied with

[TABLE]

which equivalently means that the logarithmic norm estimate is strictly negative. The interval radii of the first $39$ modes, separated above by a horizontal line, of the block $\mathcal{W}_{\text{s}}^{loc}$ are given in the eigenbasis of the finite-dimensional derivative at $z_{0}^{s,loc}$ . These coordinates are being written by the program into the file Q.in. The actual running time of the program computing this set is again negligible compared to the last step, i.e., the numerical integration. As before, we therefore do not report it.

The file stable_box_log.txt that can be found online [13] contains more detailed numerical data related with this computation.

Result of the numerical integration.

This is the third and final step of the proof of the main theorem from Section 4. Some technical parameters of the numerical integration are chosen in the following way:

•

The Taylor method order is $16$ ,

•

the constant time step is $0.002$ ,

•

the Galerkin projection dimension $m$ is $15$ , the dimension used by the FFT-algorithm $N$ is $32$ , see also the description in Section 5.

•

The whole integration process took 3m55s on an Intel®Core™i7-4610M CPU @ 3.00GHz x 4, compiled with gcc version 4.9.2 (64-bit Ubuntu 4.9.2-10ubuntu13) and with flags -O2 -frounding-math.

•

In contrast, if instead of the FFT-algorithm we used the direct algorithm, then the whole integration process took 11m55s.

We note that the overall speedup of the FFT algorithm over the direct algorithm does not seem to be overwhelming in this particular example, we guess it is due to the fact of moderate Galerkin approximation dimension ( $15$ ) used. We also guess that some parts of our code may not be optimized, which results in large constants multiplying the total cost of the FFT algorithm (17). We expect much larger speedups in case of larger Galerkin approximation dimensions. We will investigate this as a future research.

The file num_integration_log.txt that can be found online [13] contains more detailed numerical data related to this computation. Specifically, the file contains

•

the vector which is the finite-dimensional part of the propagated polynomial bounds, saved at each step of the integration,

•

each $100$ steps of the integration the full infinite-dimensional polynomial bounds,

•

at each step of the integration information at which coordinates the entry of $\mathcal{W}_{0}$ into $\mathcal{W}_{\text{s}}^{loc}$ has not yet been attained.

The infinite-dimensional bounds at time $T=3.02$ after performing $1510$ numerical integration steps are

[TABLE]

The bounds above are within the set $\mathcal{W}_{\text{s}}^{loc}$ , i.e., within the basin of attraction of the stable fixed point. This implies that the third step of the proof of the main theorem from Section 4 succeeded.

A.2 Heteroclinic to the global minimizer

We present the analogous numerical data from the computer-assisted existence proof of the heteroclinic between the homogeneous state and the global minimizer.

Unstable isolating block $\mathcal{W}_{u}\subset H^{\prime}$ around zero

This is the first step of the proof of the main theorem from Section 4). We mark with text in bold the numerical values for the eigendirection $e_{2}$ along which the unstable manifold of the homogeneous state connects to the stable global minimizer ( $\mathcal{W}_{u}$ is stretched along this direction).

[TABLE]

The infinite-dimensional cone condition (18) is satisfied with

[TABLE]

The file unstable_box_log.txt which can be found online [13] contains more detailed numerical data related with this computation.

Stable isolating block $\mathcal{W}_{\text{s}}^{glob}\subset H^{\prime}$

(The second step of the proof of the main theorem from Section 4) $\mathcal{W}_{\text{s}}^{glob}$ is centered at a numerically computed approximation to the actual global stable equilibrium loaded from the provided file fixedPointGlobal.in.

Observe that although the dimension of the finite-dimensional subtail here is larger ( $300$ ), where by finite-dimensional subtail we mean the distinguished part of the whole tail $T$ , which is represented by explicit intervals. The diameters of $\mathcal{W}_{\text{s}}^{glob}$ are now considerably larger, as well as the $\varepsilon$ cone condition value, than in the case of the local minimizer. This is due to the fact that the global minimizer is stronger attracting than the local minimizer in the sense that the maximal eigenvalue of the Jacobian matrix at the global minimizer is further away from zero than the maximal eigenvalue at the local minimizer.

[TABLE]

The infinite-dimensional cone condition (18) is satisfied with

[TABLE]

which equivalently means that the logarithmic norm estimate is negative.

The file stable_box_log.txt that can be found online [13] contains more detailed numerical data related with this computation. As in the case of the local minimizer, the interval radii of the first $39$ modes (separated above with a horizontal line) of $\mathcal{W}_{\text{s}}^{glob}$ are given in the eigenbasis of the finite-dimensional derivative at $z_{0}^{s,glob}$ .

Result of the numerical integration

This is the third and final step of the proof of the main theorem from Section 4. Some technical parameters of the numerical integration are chosen in the following way:

•

The whole integration process took 6m34s on an Intel®Core™i7-4610M CPU @ 3.00GHz x 4, compiled with gcc version 4.9.2 (64-bit Ubuntu 4.9.2-10ubuntu13) and with flags -O2 -frounding-math.

•

In contrast, if instead of the FFT-algorithm we used the direct algorithm, then the whole integration process took 17m29s.

The infinite-dimensional bounds at time $T=4.62$ after performing $2310$ numerical integration steps are

[TABLE]

The bounds above are within $\mathcal{W}_{\text{s}}^{glob}$ , i.e., the basin of attraction of the stable fixed point, which means that the third step of the proof of the main theorem from Section 4 succeeded.

Bibliography66

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] M. Bahiana and Y. Oono. Cell dynamical system approach to block copolymers. Physical Review A , 41:6763–6771, 1990.
2[2] F. A. Bartha and W. Tucker. Fixed points of a destabilized Kuramoto–Sivashinsky equation. Applied Mathematics and Computation , 266:339 – 349, 2015.
3[3] P. W. Bates and P. C. Fife. The dynamics of nucleation for the Cahn-Hilliard equation. SIAM Journal on Applied Mathematics , 53(4):990–1008, 1993.
4[4] D. Blömker, B. Gawron, and T. Wanner. Nucleation in the one-dimensional stochastic Cahn-Hilliard model. Discrete and Continuous Dynamical Systems, Series A , 27(1):25–52, 2010.
5[5] X. Cabré, E. Fontich, and R. de la Llave. T. parameterization method for invariant manifolds i: manifolds associated to non-resonant subspaces. Indiana Univ. Math. J. , 52:283–328, 2003.
6[6] J. W. Cahn and J. E. Hilliard. Free energy of a nonuniform system I. Interfacial free energy. Journal of Chemical Physics , 28:258–267, 1958.
7[7] R. Castelli, M. Gameiro, and J.-P. Lessard. Rigorous numerics for ill-posed PD Es: periodic orbits in the boussinesq equation. Preprint.
8[8] R. Choksi, M. Maras, and J. F. Williams. 2D phase diagram for minimizers of a Cahn-Hilliard functional with long-range interactions. SIAM Journal on Applied Dynamical Systems , 10(4):1344–1362, 2011.

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Code & Models

Videos

Keywords:

AMS classification:

1 Introduction

1.1 The diblock copolymer model description and the main result

Theorem 1.1** (Existence of heteroclinic connections).**

1.2 Framework for proving constructively connecting orbits in parabolic PDEs

1.3 Structure of the paper

2 Problem formulation and basic setting

Definition 2.1** (Sequence space with algebraic coefficient decay).**

Definition 2.2** (Galerkin approximation).**

3 Isolating blocks, self-consistent bounds, and the cone condition

3.1 Isolating blocks and cone conditions in finite dimensions

Definition 3.1** (hhh-set [65, Definition 1]).**

Definition 3.2** (Horizontal disk [58, Definition 10]).**

Definition 3.3** (hhh-set with cones [65, Definition 8]).**

Definition 3.4** (Cone condition [65, Definition 10]).**

Definition 3.5** (Isolating block).**

Lemma 3.6** (Graph representation [62, Lemma 5]).**

Definition 3.7** (Hyperbolic fixed point).**

Theorem 3.8** (Invariant manifolds as disks [62, Theorem 26]).**

3.2 Self-consistent bounds

Definition 3.9** (Self-consistent bounds).**

3.3 Cone conditions in infinite dimensions

Definition 3.10** (Self-consistent bounds and hhh-sets).**

Theorem 3.11** (Consequences of the cone condition).**

Lemma 3.12** (Verifying the cone condition).**

4 Proof of the main result

Theorem 4.1** (Rescaled version of the main theorem).**

Proof.

5 Rigorous integration algorithm

Definition 5.1** (Normalized derivative).**

Remark 5.2** (Recursion for normalized derivatives).**

Remark 5.3** (Computational cost of evaluating the convolution).**

Remark 5.4** (Computational cost of the convolution via the fast Fourier transform).**

Remark 5.5**.**

6 Verifying the cone condition in infinite dimensions

7 Software

8 Conclusion and future work

9 Acknowledgments

Appendix A Numerical data from the computer proof

A.1 Heteroclinic to the local minimizer

Unstable isolating block Wu⊂H′\mathcal{W}_{u}\subset H^{\prime}Wu​⊂H′

Stable isolating block Wsloc⊂H′\mathcal{W}_{\text{s}}^{loc}\subset H^{\prime}Wsloc​⊂H′

Result of the numerical integration.

A.2 Heteroclinic to the global minimizer

Unstable isolating block Wu⊂H′\mathcal{W}_{u}\subset H^{\prime}Wu​⊂H′ around zero

Stable isolating block Wsglob⊂H′\mathcal{W}_{\text{s}}^{glob}\subset H^{\prime}Wsglob​⊂H′

Result of the numerical integration

Theorem 1.1 (Existence of heteroclinic connections).

Definition 2.1 (Sequence space with algebraic coefficient decay).

Definition 2.2 (Galerkin approximation).

Definition 3.1 ( $h$ -set [65, Definition 1]).

Definition 3.2 (Horizontal disk [58, Definition 10]).

Definition 3.3 ( $h$ -set with cones [65, Definition 8]).

Definition 3.4 (Cone condition [65, Definition 10]).

Definition 3.5 (Isolating block).

Lemma 3.6 (Graph representation [62, Lemma 5]).

Definition 3.7 (Hyperbolic fixed point).

Theorem 3.8 (Invariant manifolds as disks [62, Theorem 26]).

Definition 3.9 (Self-consistent bounds).

Definition 3.10 (Self-consistent bounds and $h$ -sets).

Theorem 3.11 (Consequences of the cone condition).

Lemma 3.12 (Verifying the cone condition).

Theorem 4.1 (Rescaled version of the main theorem).

Definition 5.1 (Normalized derivative).

Remark 5.2 (Recursion for normalized derivatives).

Remark 5.3 (Computational cost of evaluating the convolution).

Remark 5.4 (Computational cost of the convolution via the fast Fourier transform).

Remark 5.5.

Unstable isolating block $\mathcal{W}_{u}\subset H^{\prime}$

Stable isolating block $\mathcal{W}_{\text{s}}^{loc}\subset H^{\prime}$

Unstable isolating block $\mathcal{W}_{u}\subset H^{\prime}$ around zero

Stable isolating block $\mathcal{W}_{\text{s}}^{glob}\subset H^{\prime}$