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

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 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.

## Authors

• 7 publications
• 1 publication
• ### Rigorous validation of a Hopf bifurcation in the Kuramoto-Sivashinsky PDE

We use computer-assisted proof techniques to prove that a branch of non-...
09/28/2020 ∙ by Jan Bouwe van den Berg, et al. ∙ 0

• ### Analysis of boundary effects on PDE-based sampling of Whittle-Matérn random fields

We consider the generation of samples of a mean-zero Gaussian random fie...
09/20/2018 ∙ by Ustim Khristenko, et al. ∙ 0

• ### Solution of the Hyperbolic Partial Differential Equation on Graphs and Digital Spaces: a Klein Bottle a Projective Plane and a 4D Sphere

In many cases, analytic solutions of partial differential equations may ...
05/01/2017 ∙ by Alexander V. Evako, et al. ∙ 0

• ### Equilibria of an anisotropic nonlocal interaction equation: Analysis and numerics

In this paper, we study the equilibria of an anisotropic, nonlocal aggre...
12/19/2019 ∙ by José A. Carrillo, et al. ∙ 0

• ### An efficient and convergent finite element scheme for Cahn--Hilliard equations with dynamic boundary conditions

The Cahn--Hilliard equation is a widely used model that describes amongs...
08/14/2019 ∙ by Stefan Metzger, et al. ∙ 0

• ### Computer-assisted proof for the stationary solution existence of the Navier–Stokes equation over 3D domains

The verified computing with assistance of computers has proved to be a p...
01/11/2021 ∙ by Xuefeng Liu, et al. ∙ 0

• ### An Entropy Equation for Energy

This paper describes an entropy equation, but one that should be used fo...
07/07/2020 ∙ by Kieran Greer, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 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 -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 . 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

 Eϵ,σ[u]=∫Ω(ϵ22|∇u|2+F(u))dx+σ2∫Ω∣∣(−Δ)−1/2(u(x)−μ)∣∣2dx. (1)

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

 ut = −Δ(ϵ2Δu+f(u))−σ(u−μ) in Ω, μ=1|Ω|∫Ωu(x)dx, and ∂u∂ν=∂Δu∂ν=0 on ∂Ω,

where the nonlinearity is defined as , and which in our situation is given by . This evolution equation is based on the gradient in the -topology and uses homogeneous Neumann boundary conditions for both  and . The diblock copolymer model (1.1) is an extension of the celebrated Cahn-Hilliard equation [6], which corresponds to the special case 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 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 . 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 , , , and , from top left to bottom right, and for total mass and domain . Each diagram indicates the structure of the equilibrium set as a function of the new parameter .

The interesting regime is the limit , and by switching to  the diagrams become easier to visualize. The horizontal straight lines in all images represent the constant function , which is clearly an equilibrium, and which we call the trivial solution. As  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 -norm of the stationary state, and the bifurcation pairs are related by a norm-preserving symmetry.

For the Cahn-Hilliard case , i.e., in the top left bifurcation diagram, pairs of solutions are created at the bifurcation points  for , and the solutions on the -th branch look qualitatively similar to the function . 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 . 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 . 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  and , 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 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 . 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 , even less is known. There are mathematical results on specific types of equilibrium solutions for two- and three-dimensional domains , 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 , for interaction lengths and , and for total mass . Then there exist heteroclinic connections between the unstable homogeneous stationary state 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 . 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 , 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 , whereas the dominant harmonic of the global minimizer is (the term of the solution’s expansion in Fourier series dominates). Also, the space of sequences with for 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

, whereas to obtain the connection to the local minimizer we stretch the block in the direction of the third eigenfunction .

### 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 . 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 . The length  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

 ut=−(uxx+λf(u))xx−λσ(u−μ), (3)

which is equivalent to (1.1) up to a rescaling of time. In this formulation, we use the parameter , and consider the partial differential equation subject to the mass and boundary constraints given by

 μ=1L∫L0u(t,x)dx and ux(t,x)=uxxx(t,x)=0 for t>0,x∈{0,L}.

As before, the nonlinearity is given by , and we consider the case 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

 u(t,x)=a0(t)+2∞∑k=1ak(t)cosπkxL, (4)

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

 a0(t)=0 for all t≥0.

Notice also that the domain length  enters the expansion. We will see later that by choosing  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

 ek(x):=cosπkxL, where k≥0,

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

 u3 = (2∑k≥1akek(x))3=8∑j1≥1∑j2≥1∑j3≥1aj1aj2aj3⋅cosπj1xL⋅cosπj2xL⋅cosπj3xL = 2∑j1,j2,j3≥1aj1aj2aj3(cosπL(j1−j2−j3)x+cosπL(j1+j2−j3)x +cosπL(j1−j2+j3)x+cosπL(j1+j2+j3)x).

The last expression can be further simplified using the fact that , as long as we recall and define  for . By reordering one further obtains

 u3=∑j1,j2,j3∈Zj1+j2+j3=0aj1aj2aj3+2∑k∈N⎛⎜ ⎜ ⎜⎝∑j1,j2,j3∈Zj1+j2+j3=kaj1aj2aj3⎞⎟ ⎟ ⎟⎠cosπkxL

If we finally substitute this expression for , together with the original series expansion for , into the differential equation (3), then one obtains after division by  and extraction of the coefficients of the basis functions  the identities

 dakdt=(−k4π4L4+λk2π2L2−λσ)ak−λk2π2L2∑j1,j2,j3∈Zj1+j2+j3=kaj1aj2aj3 for% all k>0. (5)

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

satisfies

 a0=0, as well as ak=a−k∈R % for all k∈Z.

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  denote the space , i.e., elements  are sequences such that . In addition, let denote the subspace of  which is defined by

 ˜H:={{ak}k∈Z∈H: % there exists a constant C≥0 such that |ak|≤C∤k∤6 for k∈Z},

where  for and . Finally, let the space  be given by

 H′:=˜H∩{{ak}k∈Z:a0=0 and ak=a−k for all k∈Z}.

Notice that all three spaces are Hilbert spaces.

From now on, all computations we perform in infinite dimensions will be constrained to the space . 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 be arbitrary. Then the -th Galerkin approximation of (5) is defined as

 dakdt=(−k2(πL)2(k2(πL)2−λ)−λσ)ak−λk2(πL)2∑|j1|,|j2|,|j3|≤nj1+j2+j3=kaj1aj2aj3, (6)

for any coefficient index . In other words, the -th Galerkin approximation is an -dimensional system of ordinary differential equations, which depends only on the coefficients .

## 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 -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 -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 -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

be arbitrary. Then for any vector

we let  denote a norm on , 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  for . Then we define the interval set  by

 [x]=n∏k=1[x−k,x+k]⊂Rn, as long as x−k≤x+k for k=1,…,n.

For any point and arbitrary , we let denote the open ball of radius  centered at , and we use the abbreviation for the unit ball centered at the origin. For a linear map , we let  denote the spectrum of . When considering product spaces such as , we write elements in the form , where  and  denote the first and second component vectors of , respectively, and the superscripts , and refer to the stable and unstable spaces respectively. Moreover, if denotes a continuously differentiable function, and if is a compact set, then we define

 [df(Z)]={M=(Mij)ni,j=1∈Rn×n∣∣ ∣∣Mij∈[infz∈Z∂fi∂xj(z),supz∈Z∂fi∂xj(z)]}.

Consider now an ordinary differential equation of the form

 z′=f(z),z∈Rn,f∈C2(Rn,Rn). (7)

Then we denote by  the solution of (7) which satisfies the initial condition . Notice that since we assume the twice continuous differentiability of , the solution of this initial value problem is in fact uniquely determined.

After these preparations, we now turn our attention to the concept of -sets. Informally speaking, an -set is a subset of Euclidean space, which becomes a cube after a suitable choice of coordinate system. In addition, the letter  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 -sets originated in the work on covering relations for multi-dimensional dynamical systems, see [65]. From the perspective of our specific application, the considered -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  appearing in the following definition of an -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 -set is a quadruple of the form which satisfies the following three conditions:

• The set  is a compact subset of ,

• the natural numbers  satisfy the identity , and

• the mapping , where we write , is a homeomorphism which satisfies

 cN(N)=¯¯¯¯¯¯¯¯¯¯¯¯¯Bu(N)×¯¯¯¯¯¯¯¯¯¯¯¯¯Bs(N).

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

 N−c:=∂Bu(N)×¯¯¯¯¯¯¯¯¯¯¯¯¯Bs(N),N+c:=¯¯¯¯¯¯¯¯¯¯¯¯¯Bu(N)×∂Bs(N),Nc:=¯¯¯¯¯¯¯¯¯¯¯¯¯Bu(N)×¯¯¯¯¯¯¯¯¯¯¯¯¯Bs(N),N−:=c−1N(N−c),N+:=c−1N(N+c),N:=c−1N(Nc).

and we define .

The above definition shows that an -set  is the product of two closed unit balls in some appropriate coordinate system. The two integers  and  are called the nominally unstable and nominally stable dimensions, respectively. The subscript  refers to the new coordinates given by the homeomorphism . Note that if we have , then one immediately obtains , and the identity implies . The next definition introduces the concept of a horizontal disk.

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

Let  denote an arbitrary -set. In addition, consider a continuous mapping , and define . Then we say that  is a horizontal disk in , if there exists a homotopy such that

 h(0,x)=bc(x) for all x∈¯¯¯¯¯¯¯¯¯¯¯¯¯Bu(N), h(1,x)=(x,0) for all x∈¯¯¯¯¯¯¯¯¯¯¯¯¯Bu(N), h(t,x)∈N−c for all t∈[0,1] and x∈∂¯¯¯¯¯¯¯¯¯¯¯¯¯Bu(N).

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

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

Let denote an arbitrary -set, and consider the associated splitting . Furthermore, let be a quadratic form which is given as the difference

 Q(x,y)=α(x)−β(y), for arbitrary (x,y)∈Ru(N)×Rs(N),

where and are positive definite quadratic forms. Then the pair  is called an -set with cones.

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

Consider an -set with cones given by , as introduced in the previous definition. In addition, let denote a horizontal disk as in Definition 3.2. Then we say that  satisfies the cone condition (with respect to ), if and only if for arbitrary points we have

 Q(bc(x)−bc(y))>0 as long as x≠y.

Recall that we have , and .

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  denote an arbitrary -set. We say that  is an isolating block for the vector field , if and only if  is a diffeomorphism, if the sets  and  are local sections for , i.e., the vector field  is transversal to , and if the following two conditions are satisfied:

• For all points there exists a constant such that for all .

• For all points there exists a constant such that for all .

In other words, the closed set  consists of all exit points for the flow associated with (7), and the closed set  is the set of entry points.

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

Let  denote an arbitrary -set with cones, and let be a horizontal disk which satisfies the cone condition as in Definition 3.4. Then there exists a Lipschitz continuous function such that .

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

Let and consider again the ordinary differential equation (7). Then  is called a hyperbolic fixed point for (7) if and if for all eigenvalues , where  is the Jacobian matrix of the vector field  at the point , and denotes the real part of the eigenvalue .

Now consider a hyperbolic fixed point  which is contained in a given set . Then we define the stable and unstable manifolds of  as

 WsZ(z0,φ) = {z∈Z∣∣∣φ(t,z)∈Z for all t≥0 and limt→+∞φ(t,z)=z0} and WuZ(z0,φ) = {z∈Z∣∣∣φ(t,z)∈Z for all t≤0 and limt→−∞φ(t,z)=z0},

respectively.

Important for our later applications will be the following theorem. It shows that for hyperbolic fixed points, one can always find an -set  with cones in such a way that the unstable manifold  is a horizontal disk in , and that the stable manifold  is a vertical disk.

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

Suppose that the point  is a hyperbolic fixed point for the ordinary differential equation (7). Moreover, let denote an open set with . Then there exists an -set  with cones such that the following hold:

• We have both and .

• The -set  is an isolating block for the vector field .

• The unstable manifold  is a horizontal disk in  which satisfies the cone condition.

• The stable manifold  is a vertical disk in  which also satisfies the cone condition.

Finally, the unstable manifold  can be represented as the graph of a Lipschitz continuous function over the unstable space of the linearization of  at , and this graph is tangent to the unstable space at . An analogous statement holds for the stable manifold .

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 , 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

 dudt=F(u), (8)

where  is defined on a suitable subspace of  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 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

 dakdt=Fk(a)=(−k4π4L4+λk2π2L2−λσ)ak−λk2π2L2∑j1,j2,j3∈Zj1+j2+j3=kaj1aj2aj3 for% all k>0, (9)

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

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 , 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  of the sequences , 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  and  denote the lower and upper bounds, respectively, for the -th coefficient , i.e., we will assume that elements of the constructed set satisfy , as long as  is sufficiently large. In detail, self-consistent bounds are defined as follows.

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

Consider a set of the form

 W⊕T={a∈H′∣∣(a1,…,am)∈W and (am+1,am+2,…)∈T}⊂H′,

which is the product of a finite-dimensional set and an infinite-dimensional tail . Then the set  forms self-consistent bounds for the evolution equation (9) if and only if the following conditions are satisfied.

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

 R∞⊃T=∏k>m[a−k,a+k].

Moreover, there exists a usually large integer  such that for all , i.e., all but finitely many of the intervals in  contain zero.

• There exists a constant  and an exponent such that  for all .

• The vector field  points inwards on the parts of the boundary of  which are associated with the components in the tail  in the sense that for all sequences we have

 Fk(x)>0 if xk=a−k, and Fk(x)<0 if xk=a+k, for all k>m.

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 . The block decomposition of the underlying Hilbert space  is given by the infinite sum

 H′=H′1⊕H′2⊕…, (10)

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

 |x|∞:=maxk∈N|xk| for all x={xk}k∈Z∈H′,

This norm was introduced in [59]. For later applications, we also define the projection  which maps  to  for all . Finally, within the set of self-consistent bounds 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 -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 establishes self-consistent bounds for the evolution equation (9). Then the set is an -set, if and only if the finite-dimensional part  is an -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 and a specific matrix , then there exists a stationary point for (9), and the unstable manifold  is a horizontal disk in which satisfies the cone condition.

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

Consider self-consistent bounds for (9), and assume that the orthogonal projection is an isolating block for the -th Galerkin approximation (6) of (9) for all . In particular, there is a finite number of unstable directions, all of them are confined to the finite dimensional part of , and all the directions corresponding to the infinite dimensional tail are stable. In addition, let  denote an infinite diagonal matrix with  for all indices  corresponding to the unstable directions and for the stable ones. Using this matrix, we assume that

 wt([dPmF(W)]TQ+Q[dPmF(W)])w>0 (11)

for all of the form , where .

Then there exists an equilibrium for (9) and the unstable manifold  is a horizontal disk in  which satisfies the cone condition, where we set . We use here an extended notion of -set as a product of a compact subset of and an infinite tail. Moreover, the manifold  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 both the Fourier series of all terms appearing in the right-hand side  of (9) and all partial derivatives in  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 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  the eigenvalues of the linear part of (9) are all negative, and therefore we have  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 and for all one can show that

 2infx∈W∣∣∣∂Fk∂xk(x)∣∣∣−∑ℓ≠ksupx∈W∣∣∣Qℓℓ∂Fℓ∂xk(x)+Qkk∂Fk∂xℓ(x)∣∣∣≥ε, (12)

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  being given by , 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 . 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  denotes a solution of the diblock copolymer model (1.1) on the one-dimensional domain , and for fixed parameters and . Furthermore, let  be a fixed interval length, and define a function  via

 v(τ,y)=u(λτL4,yL) % for all τ≥0 and y∈(0,L).

Then one can easily verify that  solves the diblock copolymer equation

 vτ=−(vyy+λL2⋅f(v))yy−σL2⋅λL2⋅(v−μ) on the domain % R+0×(0,L),

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

It turns out that for the parameter combinations given in Theorem 1.1 it is most convenient to choose the new domain length , 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  with , for interaction lengths and , and for total mass . Then there exist heteroclinic connections between the unstable homogeneous stationary state , which is indicated by a yellow square in Figure 4, and

1. Each of the two local energy minimizers which are indicated by yellow circles, as well as

2. 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 . This equilibrium is indicated by a yellow square in Figure 4. In addition, let denote one of the two stable local energy minimizers of (1), which are indicated by a yellow circle in Figure 4. Let 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

 Wlocs⊂H′

denote an infinite-dimensional stable isolating block which contains the stable steady state , and which forms self-consistent bounds. We also let

 Wglobs⊂H′

denote the infinite-dimensional stable isolating block which contains the stable steady state

As long as it also satisfies the cone condition, the existence of stable steady states in  and 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

 Wu⊂H′

denote the infinite-dimensional unstable isolating block which contains the known unstable equilibrium , 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 . 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 in such a way that  is an unstable isolating block for (9). For this, the set  is centered at the origin, which represents our known unstable equilibrium  for (9). More specific information regarding the specific choice of  is presented in Appendix A.1, which contains verified numerical data from the program. These results show that the topological structure of  is exactly of the form shown in Figure 6. In addition, the set  is verified to satisfy the cone condition (12), using the method presented in Section 6. Now Theorem 3.11 implies that the unstable manifold  is a horizontal disk in , and therefore the unstable manifold at the origin crosses one of the faces of . Let  denote such a face, see also the diagram in Figure 7. The whole construction of  relies on the fact that the unstable equilibrium  is given by the origin. We can therefore determine the eigenfunctions of the linearization of (9) at  explicitly, and this easily shows that all eigenfunctions are just the basis functions  for . Furthermore, recall from our brief discussion in the introduction that for our choice of  the two unstable eigendirections correspond to the eigenfunctions  and . In order to make the subsequent parts of the proof computationally faster, the set  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  (we denote the face of in this direction by ) for the case of the connection to the local minimizer. Additionally, the set  is chosen in such a way that it is much more elongated in the third eigenfunction direction for the case of the connection to the global minimizer (we denote the face of in this direction by ). See again Figure 7.

Second step: Similar to the first step, we use interval arithmetic to construct self-consistent bounds in such a way that  is a stable isolating block for (9). The block  is centered at a numerically computed approximation to the actual stable equilibrium , and the topological structure of  is as indicated in Figure 5. Analogously, the block is centered at a numerically computed approximation to the stable equilibrium  – 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 , respectively. As a consequence, the block decomposition (10) of  takes on the following form. The first  coordinates correspond to the eigenbasis of the -th Galerkin approximation (6), while the remaining coordinate directions are given by the standard basis functions  for . In other words, we have to take into account that the Jacobian at , 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 in , and in . In addition, we obtain that the complete sets , are contained in the basin of attraction of , and