Rigorous numerics for nonlinear heat equations in the complex plane of time

In this paper, we introduce a method for computing rigorous local inclusions of solutions of Cauchy problems for nonlinear heat equations for complex time values. The proof is constructive and provides explicit bounds for the inclusion of the solution of the Cauchy problem, which is rewritten as a zero-finding problem on a certain Banach space. Using a solution map operator, we construct a simplified Newton operator and show that it has a unique fixed point. The fixed point together with its rigorous bounds provides the local inclusion of the solution of the Cauchy problem. The local inclusion technique is then applied iteratively to compute solutions over long time intervals. This technique is used to prove the existence of a branching singularity in the nonlinear heat equation. Finally, we introduce an approach based on the Lyapunov–Perron method for calculating part of a center-stable manifold and prove that an open set of solutions of the Cauchy problem converge to zero, hence yielding the global existence of the solutions in the complex plane of time.


Introduction
In this paper, we consider solutions of a nonlinear heat equation with the following setting: x ∈ (0, 1), t>0, u(0, x) = u 0 (x), x ∈ (0, 1) (1.1) with the periodic boundary condition in x variable. Here the subscripts denote the derivatives in the respective variables and u 0 (x) is a given initial data. It is wellknown that solutions of the nonlinear heat equation (1.1) may blow up in finite time, that is, the L ∞ -norm of the solution tends to ∞ as t tends to a certain B < ∞. Such a B is called the blow-up time of (1.1). There are plenty of studies for blow-up problems of nonlinear heat equations. One can consult previous studies by, e.g., [7,21] and references therein. We extend the time variable of (1.1) into the complex plane. Changing the t variable into z, the aim of this paper is to study dynamics of the complex-valued nonlinear heat equation under the periodic boundary condition in x with initial data u(0, x) = u 0 (x). Here, the subscript z denotes the complex-derivative with respect to z and Re/Im denote the real/imaginary part of complex values, respectively. In his pioneering works [14,15] for this problem, Masuda considered the solution of (1.2) under the Neumann boundary condition. If the initial data u 0 (x) is sufficiently close to a real constant, he proved the existence of two analytic solutions, one of which is defined on the shaded domain of Fig. 1a and the other on its mirror image with respect to the real axis in Fig. 1b. Masuda also proved that the initial data is a constant if the solution agrees in the intersection of the two domains with Re(z) > z B , where z B ≡ B denotes the blow-up time (see Fig. 1b). It follows that, if u is the solution in the shaded domain and v is the one in the mirror image of the shaded domain, then v is the complex conjugate of u. Therefore, one can conclude that there exists multiple solutions if the imaginary part of u and v are non-trivial. It indicates that the blow-up time is expected to be a branch point 1 , but it is not shown that the blowup happens at an isolated point in the complex plane.
Following Masuda, Cho et al. [4] numerically tested dynamics of the solution of (1.2) under the periodic boundary condition. They have shown that the solution may converge to the zero function on a straight path θ def = {z ∈ C : z = te iθ , t ≥ 0} where θ ∈ (−π/2, π/2) and i = √ −1. They have also shown that the solution of (1.2) may have only one singularity on the real axis, which branches the Riemann surface of complex-valued solutions. As conclusions of their study, two conjectures have been proposed.  [14,15]. Here, 0 < θ < π/2 and z B denotes the blow-up time. b The shaded domain and its mirror-image about the real axis: intersection of two domains is plotted in dark gray • The analytic function defined by the nonlinear heat equation (1.2) has branching singularities and only branching singularities, unless it is constant in x, • The nonlinear Schrödinger equation, which is the case of θ = ±π/2, is globally well-posed for any real initial data, small or large.
The contribution of the current paper is to give two computer-assisted proofs for the complex-valued nonlinear heat equation (1.2). Although these results do not fully solve the above conjectures, they do provide a mathematical proof of the dynamical behavior of the solution numerically obtained by Cho et al. [4]. Additionally they also extend the results by Masuda [14,15] in terms of treating the periodic boundary conditions.
Our first result shows the existence of a branching singularity, which appears at the blow-up point Cho et al. [4] calculated to be approximately z B ≈ 0.0119. This theorem partially answers the first conjecture, that is it proves that the analytic function defined by (1.2) has at least one branching point. If we assume the second conjecture, then we would be able to conclude that there is a real blow-up point z B in the interval (0.0116,0.0145). However even then it is not clear that the set of branching singularities is a single point, and not a branch cut along an interval, cf. [24]. Our second result agrees with Masuda's work for the case of periodic boundary conditions without assumption of closeness to a constant for the initial data. Theorem 1.2 (Global existence) For θ = π/3, π/4, π/6 and π/12, setting a straight path θ : z = te iθ (t ≥ 0) in the complex plane of time, the solution of the complexvalued nonlinear heat equation (1.2) under the periodic boundary condition with the specific initial data u 0 (x) = 50(1 − cos(2π x)) exists globally in t and converges to the zero function as t → ∞.
The proofs of these results are obtained in Sect. 7 by using rigorous numerics, via a careful blend of functional analysis, semi-group theory, numerical analysis, fixed point theory, the Lyapunov-Perron method and interval arithmetic. It is worth mentioning that the study of finite time blow-up in ODEs [6,16,26] and global existence of solutions in PDEs [2,18,29] are beginning to be studied with the tools of rigorous numerics. Furthermore, several approaches of rigorous integration for evolution PDEs have been introduced in, e.g., [5,10,17,30,31]. We expect that these rigorous integrators will be applied to understand dynamics of various evolution PDEs soon. In recent works, we have applied the rigorous integrator developed in this paper to study the dynamics of (1.2) in further detail, proving the existence of heteroclinic orbits and other global solutions [8,9].
It is also important to note that the techniques of rigorous numerics presented in this paper are by no means applicable only to nonlinear heat equations with the quadratic nonlinear term. These techniques can be easily applied to more general class of nonlinearities, such as u p when p ≥ 2 is an integer. Moreover, we believe that our approach should be extended in principle to more general PDE models involving e u or even a derivative type nonlinear term ∂ x (u p ).
The present paper is organized as follows. In Sect. 2 we set up the fixed point operator to solve the Cauchy problem (1.2) on a straight path in the complex plane. Using the Fourier series, we derive an infinite-dimensional system of differential equations and define a formulation of the fixed point operator to rigorously provide an enclosure of the solution on a short time interval. Section 3 is devoted to defining a solution map operator, which provides a solution of the linearized Cauchy problem with arbitrary forcing term. Such a solution map operator is defined by what is called the evolution operator. We show how we verify the existence of such an evolution operator by using rigorous numerics and give a computable estimate for a uniform bound of the evolution operator. Using the solution map operator, we derive in Sect. 4 a sufficient condition to have a fixed point for the fixed point operator. Theorem 4.1 provides rigorous enclosure of the solution of the Cauchy problem. In Sect. 5, we introduce a method for applying iteratively the local inclusion approach to compute solutions over long time intervals. This technique is then used to prove Theorem 1.1. In Sect. 6, we introduce an approach based on the Lyapunov-Perron method in order to calculate part of a center-stable manifold, which is used to prove that an explicit open set of solutions of the given Cauchy problem converge to zero, hence yielding the global existence in time of the solutions. This proves our second main result, namely Theorem 1.2. The numerical results are presented in details in Sect. 7.

Setting-up the fixed point operator to solve the Cauchy problem
In this section, we derive the formulation of the fixed point operator which is used to provide a rigorous enclosure for the solution of the Cauchy problem on a short time interval. Taking the straight path θ = {z ∈ C : z = te iθ , t ≥ 0}, the complex-valued nonlinear heat equation (1.2) is transformed into the equation with periodic boundary conditions. Here, we set θ ∈ (−π/2, π/2) and the initial data u 0 (x) = 50(1 − cos(2π x)). We expand the unknown function into the Fourier series: where we have set ω = 2π . Plugging (2.2) in the initial-boundary value problem (2.1), we have the following infinite-dimensional system of ODEs: where " * " denotes the discrete convolution product defined by for bi-infinite sequences b = (b k ) k∈Z and c = (c k ) k∈Z , and ϕ is defined by An important and useful feature of 1 is that it is a Banach algebra under discrete convolution, namely (2.6) To determine the Fourier coefficients (a k (t)) k∈Z of the solution of the Cauchy problem (2.3), we use a numerically computed approximate solution with the maximal wave number N , say (ā k (t)) |k|≤N . Let be an approximation of a(t). We will rigorously include the Fourier coefficients in the neighborhood of numerical solution defined by Define the Laplacian operator L acting on a sequence of Fourier coefficients as The domain of the operator L is denoted by D(L) def = a = (a k ) k∈Z : k∈Z k 2 |a k | < ∞} ⊂ 1 . We define an operator acting on a∈D Then we consider the Cauchy problem (2.3) as the zero-finding problem F(a) = 0 with the initial condition a(0) = ϕ. It is obvious that a solves (2.3) if F(a) = 0 with a(0) = ϕ holds. Let us define the following operator: where A ϕ : X → D ∩ X ⊂ X is a solution map operator corresponding to the linearized problem of (2.3) around the approximate solutionā, which will be explicitly defined in the next section. As shown in Remark 3.1, this operator is alternative form of the simplified Newton operator. We expect that the simplified Newton operator (2.10) has a fixed pointã ∈ B J (ā, ) such thatã = T (ã).

The solution map operator
The simplified Newton operator (2.10) is characterized by the solution map operator. We define the solution map operator A φ using an evolution operator is the solution of the following homogeneous initial value problem (IVP): Here, s is also a parameter since the evolution operator U (t, s) is a two parameter family of bounded linear operators. We also note that ( Setting the initial sequece b(0) = φ of (3.2), this leads to the definition of the solution map operator A φ : X → D ∩ X ⊂ X given by In particular, This is a linearized problem of (2.3). Setting the initial data (e.g., b(0) = φ), we will validate that the problem has a unique solution that is continuous on J satisfying for ϕ 1 , ϕ 2 ∈ 1 , g 1 , g 2 ∈ X . Then, the simplified Newton operator defined in (2.10) satisfies for a ∈ D ∩ X with a(0) = ϕ. We note that this form a − A 0 F(a) is defined only on D ∩ X but the simplified Newton operator (2.10) can be defined on X . Moreover, from the property of the operator A ϕ : X → D ∩ X ⊂ X , which will be shown in Remark 3.6, the fixed point satisfiesã ∈ D ∩ X if such a fixed point is obtained.
More precisely, we will derive a hypothesis of existence of the evolution operator in Theorem 3.3, which can be verified via rigorous numerics. If such a hypothesis holds, the explicit value of W h > 0 is also obtained. The constant W h > 0 provides a uniform bound for the bounded linear operator norm of U (·, ·) over the simplex Here, we use the notation B( 1 ) to denote the set of all bounded linear operators on 1 with corresponding bounded operator norm · B( 1 ) .

Remark 3.2
It is important to notice that for a givenā(t) ∈ 1 and h small enough, there exists the solution of the linearized problem (3.1) (e.g. see [19]). This in turns yields the existence of the evolution operator U (t, s) for all (t, s) ∈ S h for a small enough step size h > 0. However, once a fixed h > 0 has been chosen, it is in general difficult to show the existence of the evolution operator U (t, s) for (t, s) ∈ S h . In the next section, we derive a hypothesis of existence of the evolution operator. We also give an explicit and computable formula for the constant W h .

Deriving a formulation for W h and existence of the evolution operator
We begin this section with some notation. Given a finite dimensional (Fourier) projection number m ∈ N and a vector φ = (φ k ) k∈Z ∈ 1 , define the projection π (m) : 1 → 1 as follows: (∞) . For example, with this notation, the discrete convolution in (3.1) can be expanded as As previously mentioned, given a step size h > 0, our goal in this section is to show the existence of the evolution operator U (t, s) of (3.1) by computing a constant W h satisfying (3.4). However, from a computational point of view, the formulation of the homogeneous non-autonomous linear IVP (3.1) is daunting due to the fact that for each fixed k ∈ Z the convolution term (ā (t) * b(t)) k involves all Fourier modes (b j (t)) j∈Z . In order to address this issue, we separate the equations by considering yet another homogeneous IVP with respect to the sequence This new decoupled formulation, while not being equivalent to (3.1), will be used to control the evolution operator associated to (3.1). Setting the initial sequence c(s) = φ (0 ≤ s ≤ t), let (c k (t)) |k|≤m and (c k (t)) |k|>m be the solution of the (2m +1)dimensional equation (3.5) and the infinite dimensional equation (3.6), respectively. We define the solution operator of the IVP of (3.5) by We call C (m) (t, s) and C (∞) (t, s) the evolution operators of (3.5) and (3.6), respectively. At this point, the existence of the operator C (∞) (t, s) is not guaranteed, but will be verified in Sect. 3.3. We extend the action of the operator C (m) (t, s) (resp. C (∞) (t, s)) on 1 by introducing the operatorŪ (m) (t, s) (resp.Ū (∞) (t, s)) as follows.
The proof of existence of the evolution operator U (t, s) of the original linearized problem (3.1) and the explicit bound W h satisfying (3.4) is presented in the following theorem.
Assume that C (∞) (t, s) exists and thatŪ (∞) (t, s) defined in (3.8) satisfies then the evolution operator U (t, s) exists and the following estimate holds Remark 3.4 Equations (3.11) and (3.12) above implicitly assume that the denominator is non-zero. However, if the denominator is zero, these are W ∞ = h andW ∞ = h 2 /2. It easily follows from the proof of Lemma 3.5.
After the proof of Theorem 3.3, we introduce in Sect. 3.2 a rigorous computational method based on Chebyshev series to obtain the finite dimensional evolution operator C (m) (t, s), which will directly yieldŪ (m) (t, s). We show how this helps computing W m > 0 satisfying (3.9). Section 3.3 is devoted to show the existence of C (∞) (t, s) (and hence the existence ofŪ (∞) (t, s)) and to show that the hypothesis (3.10) holds. The proof of Theorem 3.3 uses the following elementary result.
Proof First, note that from (3.10) Second, note that Third, and hence, it follows that Proof of Theorem 3.3 First note that for |k| ≤ m, the system of differential equations (3.1) is described by the non-homogeneous equation with the initial condition b k (s) = φ k for |k| ≤ m. Consider the homogeneous equation (3.5) and denote by C (m) (t, s) the solution of the variational problem. We show in Sect. 3.2 how to rigorously compute C (m) (t, s). As before, letŪ (m) (t, s) be the extension of the action of the operator C (m) (t, s) on 1 . Using the evolution opera-torŪ (m) (t, s) : 1 → 1 , we can integrate the system (3.19) to obtain the integral equation (3.20) Here, for |k| ≤ m, where the last equality holds because the index k 1 is never equal to zero. Combining (3.9) and (3.20), and using the property (2.6), it follows that Next, for the case of |k| > m, we rewrite the system of differential equations (3.1) as (3.22) with the initial condition b k (s) = φ k for |k| > m. Note that by assumption the evolution operator C (∞) (t, s) associated to (3.6) exists. DefineŪ (∞) (t, s) as in (3.8).
Using that operator, the system (3.22) can be re-written as For |k| > m, the following holds Combining (3.10), (3.23) and using the property (2.6), it follows that (3.25) By assumption (3.14), κ = 1 − 4W mW∞ ā (s) 2 X > 0 and using (3.25) yields that Note that if the hypothesis (3.14) holds, it guarantees the existence of the solution of the finite part of (3.1). Using the inequalities (3.16) and (3.17) in Lemma 3.5 and (3.26), the tail mode (3.24) is bounded by

Remark 3.6 For arbitrary forcing term
2) in the sense of classical solution (see, e.g., [19]). This implies b ∈ D ∩ X . Then, the solution map operator A has a smoothing property from X to D ∩ X .
Having derived a computable formulation in (3.15) for the constant W h > 0 which provides a uniform bound for the norm of evolution operator U (·, ·) over S h , we now turn to the question of computing the bound W m , which is required to defined W h .

Computation of the bound W m
The goal of this section is to present a rigorous computational approach to obtain the constant W m satisfying (3.9), that is a uniform bound for the operator norm ofŪ (m) over the simplex S h . From (3.7), definingŪ (m) (t, s) boils down to the computation of C (m) (t, s), that is the evolution operators of the (2m + 1)-dimensional equation (3.5), which can be written component-wise as for 0 ≤ s ≤ t ≤ h, and where δ k, j denotes the Kronecker's delta. Denote by A(t) the matrix representation of the right-hand side of (3.28) acting on the coefficients c k, j (t, s) |k|,| j|≤m .
The evolution matrix C (m) (t, s) can be decomposed as where (t) ∈ M 2m+1 (C) is the principal fundamental solution and solves and where (s) = (s) −1 is the solution of the adjoint problem Having introduce the set-up, the computation of W m is twofold. First, using the approach of [12], we use the tools of rigorous numerics and Chebyshev series expansion to compute the fundamental matrix solutions (t) and (t) satisfying (3.30) and (3.31), respectively. Second, using interval arithmetic and the Chebyshev series representations of (t) and (t), we obtain W m such that where · 1 represents the matrix norm induced by the 1 norm on C 2m+1 . By construction, the constant W m obtained computationally in (3.32) satisfies (3.9). The remainder of this section is devoted to the computational method for obtaining the matrix solutions (t) and (t).

Rigorously computing 8(t) and 9(t) via Chebyshev series
The Chebyshev polynomials T k : has a unique representation as an absolutely and uniformly convergent series v(t) = ∞ k=0 a k T k (t) (e.g. see [27]). The more regular the function v is, the faster the decay rate of its Chebyshev coefficients.
Fixing a Fourier projection number N and a Chebyshev projection number n, let us assume that we have numerically computed an approximate solution of the initialboundary value problem (2.1) in the form with ω = 2π . In practice, we perform this task by considering a Galerkin approximation of (2.3) of size 2N + 1 and using MATLAB software system Chebfun (e.g. see [20]) Assume that the initial condition of (3.34) is given by c k (−1) = b k . In practice, we will solve rigorously 2m + 1 IVPs with initial conditions (b k ) m k=−m = e j ∈ C 2m+1 the canonical basis vectors. Rewriting the system (3.34) as an integral equation results in For each k, we expand c k (t) using a Chebyshev series, that is (3.38)

ψ ,k (c)T (s) ds
and this results (e.g. see in [12]) in solving Hence, for > 0 and |k| ≤ m, we aim at solving Finally, the problem that we solve is Define the operators (acting on Chebyshev sequences) by Using the operators T and , we may write more densely for the cases > 0 and |k| ≤ m Hence, , denote the Chebyshev- The Banach space in which we prove the existence of the solutions of f = 0 is given by The following result is useful to perform the nonlinear analysis when solving f = 0 in X m ν . We omit the elementary proof which can be mimicked from the one of Lemma 3 in [28].
where r is a radius to be determined. Let us now define the operator A. Given n, a finite number of Chebyshev coefficients used for the computation of c, denote by f (n,m) : C n(2m+1) → C n(2m+1) the finite dimensional projection used to computec ∈ C n(2m+1) , that is First consider A † an approximation for the Fréchet derivative D f (c): The following Newton-Kantorovich type theorem (for linear problems posed on Banach spaces) is useful to show the existence of zeros of f .
Proof We omit the details of this standard proof. Denote κ

and by construction of the operators A and
A † , it can be shown that A is an injective operator. By injectivity of A, we conclude that there exists a uniquec ∈ B r (c) such that f (c) = 0.
Given a Fourier projection dimension m, we apply Theorem 3.9 to compute the solution of 2m + 1 problems of the form f = 0 given in (3.42) with initial conditions Finally we can define the fundamental matrix solution (t) ∈ M 2m+1 (C) satisfying (3.30) as Using a similar construction, we can construct (s) ∈ M 2m+1 (C) satisfying (3.31). The rest of this section is dedicated to the explicit construction of the bounds Y 0 , Z 0 and Z 1 of Theorem 3.9.
The bound Y 0 Recalling the definition of the quadratic term N k of f in (3.39), given the This implies that the term f (c) = ( f j (c)) j ∈I has only finitely many nonzero terms.
The bound Z 0 The computation of the bound Z 0 satisfying (3.44) requires defining the operator which action is given by Using interval arithmetic, compute Z 0 such that The bound Z 1 . For any c ∈ B 1 (0), let which is given component-wise by Next, for the cases 0 ≤ < n and |k| ≤ m, we These bounds will then be used to define Z 1 . First, given c ∈ B 1 (0) and |k| ≤ m, The next step is to obtainẑ ,k for = 1, . . . , n −1. This task involves understandinḡ a * c (∞,m) whose components can interpreted as For 0 < < n and |k| ≤ m, the term (ā * c (∞,m) ) ,k ∈ C can then be thought of as a linear functional acting on c = (c j ) j ∈I ∈ X m ν . Using that representation, we get that for all c ∈ X m ν with c X m ν ≤ 1, which can easily be computed using interval arithmetic. For the cases 0 ≤ < n and |k| ≤ m, this leads to the bound for all c ∈ B 1 (0), where |T | denotes the operator with component-wise absolute values. We are ready to obtain the bound Z 1 . Given c ∈ B 1 (0), we use Lemma 3.8 to conclude that Hence, by construction satisfies (3.45).

Generation of the evolution operatorŪ (∞) (t,s) on 1
Finally, in this section, we verify the existence of the evolution operator C (∞) (t, s) of the infinite dimensional equation (3.6). This in turns will verify the existence of the evolution operatorŪ (∞) (t, s) generated on 1 . Moreover, we derive an estimate ofŪ (∞) (t, s) defined in (3.8) using the bounded operator norm on 1 , which is the hypothesis (3.10) of Theorem 3.3.
where the domain of the operator L is D(L) = a ∈ 1 ∞ : La ∈ 1 ∞ in this case. Denoting L = e iθ L, it is easy to see that L : D(L) ⊂ 1 ∞ → 1 ∞ is a densely defined closed operator on 1 ∞ . The operator L then generates the semigroup on 1 ∞ (e.g. see [19]), which is denoted by e Lt t≥0 . Furthermore, the action of such a semigroup e Lt t≥0 can be naturally extended to 1 by In the following, unless otherwise noted, we consider the semigroup e Lt t≥0 on 1 as defined in (3.52). We also have the following estimate: (3.53) We re-write (3.6) as the Cauchy problem in 1 for any initial sequence c (∞) (s) = φ (∞) . Showing the existence of the solution of (3.54), the proof of existence of the evolution operatorŪ (∞) (t, s) is given by the following theorem.

Theorem 3.10
For the infinite-dimensional system of differential equations (3.54) with a fixed initial time s, there exists a unique solution that solves the integral equation in Furthermore, the evolution operatorŪ (∞) (t, s) exists and the following estimate holds Proof For a fixed s ≥ 0, let us define a map P : X → X acting on the c (∞) (t) as and let a function space X ∞ be defined by We prove that the map P becomes a contraction mapping under the distance d on ∈ X ∞ , we have using (3.53) and the property (2.6) Since β > 1, P becomes a contraction mapping on X ∞ . This yields that the solution of (3.54) uniquely exists in X ∞ , which satisfies 1 , it follows from (3.55) using (3.53) and the property (2.6), that Grönwall's inequality yields Then, we conclude that the following inequality holds for any φ (∞) ∈ 1 , where W (∞) (t, s) is defined in (3.10).

Rigorous construction of the solution map operator
Recall that the definition of the solution map operator A φ in (3.3) requires proving the existence of the evolution operator U (t, s), which is done by verifying the main three hypotheses of Theorem 3.3, namely (3.9), (3.10) and (3.14). Hypothesis

Local inclusion in a time interval
In this section, we present a sufficient condition whether a fixed point of the simplified Newton operator (2.10) exists and is unique in B J (ā, ) defined by (2.7). Such a sufficient condition can be rigorously checked by numerical computations based on interval arithmetic.
Assume also thatā ∈ X and any a ∈ B J (ā, ) satisfies The proof of this theorem is based on Banach fixed point theorem.
Proof On the basis of Banach fixed point theorem, we prove that the simplified Newton operator T defined by (2.10) becomes a contraction mapping on B J (ā, ). It is sufficient to show that the following two conditions hold: Firstly, for a sequence a ∈ B J (ā, ) , we have using (2.10) Let z def = a −ā. The first term of (4.2) follows Thus, (4.2) is represented by where A z(0) is the solution map operator defined in Sect. 3 and Taking 1 -norm of g, we have using the property (2.6) where δ satisfies sup s∈J (F(ā)) (s) 1 ≤ δ. We remark that the last inequality in (4.6) is obviously overestimated, but we estimate it in this way to simplify the later notations in the second part of this proof. Since a ∈ B J (ā, ), z X ≤ holds. Finally, (4.4) is bounded by using the uniform bound W h (3.4) discussed in the previous section as sup t∈J k∈Z where ε is the upper bound of the initial error such that z(0) 1 ≤ ε. From the assumption f ε ( ) ≤ , T (a) ∈ B J (ā, ) holds for any a ∈ B J (ā, ). Secondly, we will show the contraction property of T . For sequences a 1 , a 2 ∈ B J (ā, ), we define the distance in B J (ā, ) as The analogous discussion above yields Then, from (4.7) and (4.8), the distance is estimated by a 2 ).
Taking κ = 2W h h , it follows κ < f ε ( )/ ≤ 1 from the assumption of theorem. It is proved that the simplified Newton operator T becomes the contraction mapping on B J (ā, ).

Remark 4.2
Our task in the practical implementation is to rigorously compute the minimum values such that f ε ( ) ≤ by using interval arithmetic.

The bound "
We show how we get the ε bound such that ϕ −ā(0) 1 ≤ ε. From (3.33) where we used the fact T (0) = (−1) . Then, using interval arithmetic, ε is given by We remark that ϕ k = 0 for |k| > N from the definition (2.4). If these are non-zero, ε should include the tail term |k|>N |ϕ k |.

The bound ı
We also show how we get the defect bound of F at the approximate solutionā. From (4.3), we recall Here, we suppose that the first derivative ofā is expressed by ,k ∈ C can be computed by an recursive algorithm (see, e.g., [13,page 34]). For |k| ≤ N , we have where w denotes a weight for the Chebyshev coefficients such that It then follows Furthermore, for |k| > N , the tail part is given by Finally, the defect bound δ is given by This seems to be an infinite sum, but thanks to the finite number of nonzero elements inā ,k , it becomes a finite sum. Then the defect bound can be rigorously computed by interval arithmetic. For the rigorous computing of the discrete convolution, one can consult a FFT based algorithm by, e.g., [11].

Time stepping scheme
Once the rigorous inclusion of Fourier coefficients is obtained, we consider extending the time interval, say time step, in which the existence of the solution is verified. For this purpose the initial sequence in the next time step is replaced by a sequence at the endpoint of the current time step (e.g., a(h) for the first time step). Replacing J = (h, 2h), we apply Theorem 4.1 for the initial-boundary value problem on the next time step and recursively repeat this process several times. We introduce such a time stepping scheme in this section.
For K ∈ N, let 0 = t 0 < t 1 < · · · < t K < ∞ be grid points of the time variable. We call J i = (t i−1 , t i ) the i th time step and let h i = t i − t i−1 (i = 1, 2, . . . , K ) be the step size. Now we assume that the solution a(t) = (a k (t)) k∈Z of (2.3) is rigorously included until J K , i.e., holds for some i > 0. In the following, we derive the error bound ε i (i = 1, 2, . . . , K ) at each endpoint of time interval, which satisfies a(t i ) −ā(t i ) 1 ≤ ε i . We call such an error estimate the point-wise error estimate.
Let us show the case of first time step. Recall that z(t) = a(t) −ā(t) (t ∈ J 1 ). At the endpoint of the time step t = t 1 , it follows from (4.4) where g(s) is defined in (4.5) and g(s) 1 ≤ 2 2 1 + δ 1 holds from (4.6). Here, the positive constant δ 1 satisfies sup s∈J 1 (F(ā)) (s) 1 ≤ δ 1 , which is given in Sect. 4.2. We assume that the initial error is bounded by z(t 0 ) 1 ≤ ε 0 , which is the same as that in the hypothesis of Theorem 4.1 in Sect. 4. Then, to derive the bound ε 1 , we need two positive constants W J 1 > 0 and W t 1 > 0 such that sup s∈J 1 U (t 1 , s) We reconsider the linearized problem (3.1) and represent the solution of (3.1) as b s (t) ≡ b(t). As the analogous discussion in Sect. 3, we split the solution b s (t) into the finite mode b s . The finite mode is given by plugging Taking the supremum norm with respect to s, it follows Next, the tail mode is given by plugging t = t 1 in (3.23) Using the bound (3.56), the property (2.6) and (3.26), we also have Taking the supremum norm with respect to s, it follows Moreover, let us define Summing up (5.2) and (5.3), we get the W J 1 bound sup s∈J 1 where

Setting s = t 1 , (5.2) and (5.3) are changed by
respectively. The analogous discussion yields the W t 1 bound by Consequently, the point-wise error is estimated by For the next time step, updating ε ≡ ε 1 , we apply Theorem 4.1 on J 2 and derive the point-wise error ε 2 at the endpoint recursively. By repeating the time stepping, we 1, 2, . . . , K ).

Remark 5.1
The point-wise error ε i can be smaller than the previous one, i.e., ε i ≤ ε i−1 holds when W t i < 1 with sufficiently small i and δ i . Furthermore, we note that the tiny numerical error may include at the end point. That is, the value of approximate solution at t = t 1 in J 1 , sayā J 1 (t 1 ), and that in J 2 , sayā J 2 (t 1 ), may be different. It is because the Chebyshev polynomial approximates a function globally in each time interval. In such a case, we should take care of the numerical error by using the following form: Hence, we should add 1 in the point-wise error, i.e., ε 1 = W t 1 ε 0 +W J 1 h 1 (2 1 + δ 1 )+ 1 .

Global existence in time
After several time steppings, we prove global existence of the solution. Our approach is based on a calculation of a center-stable manifold arising from the Cauchy problem (2.3).

Calculating part of a center-stable manifold
Starting from the PDE (2.1) we consider the system of differential equations in 1 given by (2.3) In this section, we use the Lyapunov-Perron method for computing a foliation of portion of the center-stable manifold of the equilibrium at a ≡ 0. A good reference for this method in ODEs is [3], and for PDEs see [23]. In [29] this method is applied to give computer-assisted proofs of the stable manifold theorem in the Swift-Hohenberg PDE. Let us define subspaces We may then rewrite (2.3) as the following system: where for a = (x c , x s ) we define: Using the convolution product, the nonlinearities may be written as: To abbreviate, let us define μ = ω 2 cos θ . One may observe that the center subspace X c is in fact invariant with respect to the nonlinear dynamics; X c is a center manifold of the zero equilibrium. We can solve (6.1) restricted to this subspace:ẋ c = e iθ x 2 c , (6.4) using the fact that the differential equation is separable. For an initial condition φ ∈ C then the solution of (6.4) is given by: Most solutions limit to zero in forward and backwards time however if Re(φe iθ ) > 0 and Im(φe iθ ) = 0, then the solution (t, φ) blows up in finite time.
For r c , r s > 0 let us define the following sets Note that if φ ∈ B c (r c ) then |φ| 1+|φ|t ≤ | (t, φ)| ≤ |φ| √ 1+|φ| 2 t 2 for t ≥ 0. Moreover | (t, φ)| ≤ r c and Re(e iθ (t, φ)) ≤ 0 for t ≥ 0, hence B c (r c ) is a forward invariant set (see Fig. 2). We aim to locally characterize the center-stable manifold at 0 as a Lipschitz graph over B c . To that end we make the following definition.
Since (t, φ) limits to zero, such a fixed point α = [α] gives us a foliation of the center stable manifold over B c (r ). By Corollary 6.4 we obtain an explicit neighborhood within which all points limit to the zero equilibrium. To prove the existence of a fixed point to our Lyapunov-Perron operator, and obtain explicit bounds, we prove the following theorem. Theorem 6.2 Fix r c , r s , ρ > 0 and define the following constants If the following inequalities are satisfied The proof is given in Sect. 6.2. , ξ ), ξ ). Then lim t→∞ X (t) = 0 and moreover

Corollary 6.3 Suppose the hypothesis of Theorem 6.2 is satisfied andα
The proof largely follows from Proposition 6.8.
where ∂ B c denotes the boundary of B c . Suppose the hypothesis of Theorem 6.2 is satisfied and there exists some We wish to show there exists some φ ∈ B c such that α(φ, ξ ) = φ 0 . Let us define To show U is contained in the α-skew image of B c × B s , it suffices to show that τ 0 = 1.

continuous curve, and we may define a function
: By way of contradiction, let us suppose that τ 0 < 1. If φ 1 ∈ int(B c ), then H maps  an open set about (φ 1 , τ 0 ) onto an open set about (φ 0 , τ 0 ). Hence the initial choice of τ 0 as the supremum above was incorrect, and τ 0 could have been made even larger. Otherwise suppose φ 1 ∈ ∂ B c . But this is contradicted by our definition of U , as seen by the following calculation: That is to say a(0) is in the center stable manifold of the equilibrium 0; lim t→∞ a(t) = 0. Remark 6.5 By isolating the ρ terms in the inequality δ 3 μ−δ 2 < ρ, we obtain the following: Hence, we need at a minimum μ > 4r c + 2r s in order to satisfy the hypotheses of Theorem 6.2. Using the quadratic formula, we obtain a necessary inequality, We can use this to define a nearly optimal ρ in terms of r c , r s and μ.

Remark 6.6
For a given θ , there is something to be said about how to select r c and r s so that the hypothesis of Theorem 6.2 is satisfied. For the best bounds, ρ should be taken as small as possible, and nearly optimal bounds are explicitly given in Remark 6.5. Furthermore, we necessarily need to take r c < 1 4 ω 2 cos θ . Below are some constants which will satisfy Theorem 6.2 and try to maximize r c .

Proof of Theorem 6.2
The proof of Theorem 6.2 is organized as follows. First in Propositions 6.8-6.10 we derive bounds on the norm of solutions to (6.7) and their dependence on initial conditions. In Proposition 6.12, we show that : B → B is a well defined endomorphism. In Proposition 6.14 we introduce a norm · Ẽ onB. In Proposition 6.15, we show that if two maps α, β ∈ B are such that α − β Ẽ is small, then their solutions to (6.7), taken with the same initial conditions but different maps α and β, will also remain close. Finally, we finish the proof of Theorem 6.2 by bounding the Lipschitz constant of using this new norm. If this Lipschitz constant is less than 1, then is a contraction mapping, and hence it will have a unique fixed point. Throughout, we will often make use of the following Grönwall type lemma below [29].
Proof Let us abbreviate x s = x s (t) = x(t, φ, ξ, α) and let us write α(φ, ξ ) = φ + α(φ, ξ ) forα ∈B. By variation of constants we have Since | (t, φ)| ≤ r c for t ≥ 0, then by (6.5) we obtain : , φ), 0))| ≤ r c + ρ x s . (6.11) Let us now assume that ξ < r s , by which there exists some positive T def = sup{t > 0 : x s (t) ≤ r s }. Hence |α( (t, φ), x s )| ≤ r c + ρr s for t ∈ (0, T ). Using the estimate in (6.3) we obtain N s (α( (t, φ) for all t ∈ [0, T ). Returning to our variation of constants formula, we have: From Grönwall's inequality, we obtain x s (t) ≤ e −(μ−δ 1 )t ξ for all t ∈ (0, T ). From this we may conclude that in fact T = +∞. By continuity the inequality also extends to the case ξ = r s . Proposition 6.9 Fix α ∈ B, φ ∈ B c (r c ), and ξ, ζ ∈ B s (r s ), and define Proof Let us abbreviate x s = x s (t) = x(t, φ, ξ, α) and z s = z s (t) = x(t, φ, ζ, α). Additionally, let us abbreviate x c = α( (t, φ), x s (t)) and z c = α( (t, φ), z s (t)). By variation of constants we have Using the notational abbreviation x 2 s = x s * x s for x s ∈ X s , we bound the norm of the difference in the integrand below By Proposition 6.8 then x s , z s ≤ r s for all t ≥ 0. So from (6.11) we obtain |x c | ≤ r c + ρr s . By the Banach algebra property of * we have and thereby we have x 2 s − z 2 s 1 ≤ 2r s x s − z s 1 for all t ≥ 0. Also note that by (6.5) then Plugging all of this into (6.12) we obtain Returning to our variation of constants formula, we have: The proposition follows from Grönwall's inequality.
To prove the proposition we will show that there exist 0 , 1 > 0, given explicitly in (6.24) and (6.16) respectively, such that Thus the difference in (6.17) can be expressed by the integral below, which we estimate in three parts: Part I By the triangle inequality, the integrand in Part I can be bounded as follows Since (t, φ) ∈ B c for all t ≥ 0 and φ ∈ B c , then | (t, φ 1 )| ≤ r c for all t ≥ 0. To continue bounding the first summand we apply inequalities (6.5) and (6.6) as below: To bound these terms, we use Proposition 6.8, inequality (6.14), and Proposition 6.10 as below: . (6.19) From the estimate in (6.14) we have We finish bounding the second summand in Part I using (6.5) and Proposition 6.8 as below Combining all of these estimates into inequality (6.18), we obtain a bound for the integrand of Part I Integrating we obtain our bound on Part I,

(6.21)
Part II To estimate the integrand in Part II, we begin by writing: In equations (6.19) and (6.20) we calculated similar bounds for these terms, by which we obtain Integrating, we obtain a bound for Part II.
Part III Using the Banach algebra property of * , the integrand in Part III may be re-expressed as follows: Using the estimates from Propositions 6.8 and 6.10 we obtain the following bound Integrating, we obtain a bound for Part III.
To show that is a contraction mapping, we will use a somewhat weaker norm.

Definition 6.13 Define the vector space below
This can be seen to be a Banach space by the same argument as in [3]. Moreover E). Note that all functions inB ρ,κ have a uniform bound on their Lipschitz constant. Hence the sequence is uniformly equicontinuous, and converges uniformly to a continuous functionα ∈ C 0 (B c × B s , X c ). Moreover, asα n (φ, 0) = 0 for all φ ∈ B c and n ∈ N, thenα(φ, 0) = 0. Furthermore, since there is a fixed bound on the Lipschitz constant for the entire sequence, thenα is Lipschitz continuous. By the same argumentα also satisfies the Lipschitz bounds in (6.5) and (6.6), henceα ∈B ρ,κ .
Additionally, let us abbreviate x c = α( (t, φ), x s (t)) and y c = β( (t, φ), y s (t)). By variations of constants we have the following: In the same manner as in (6.12), we bound the difference of the nonlinearities: where we used Proposition 6.8 to obtain x s , y s ≤ r s , and (6.11) to obtain |x c | ≤ r c + ρr s . To bound |x c − y c | we calculate further, using (6.5) and the definition of · Ẽ to obtain the following: Combining these estimates we obtain: We may continue calculating using variation of constants: By Lemma 6.7, we have that: The proposition follows.
We can now show that is a contraction map.
Proof of Theorem 6.2 By the Proposition 6.12 there exists a κ > 0 such that : B ρ,κ → B ρ,κ is an endomorphism. Define the function i c : B c × B s → X c by i c (φ, ξ ) = φ, and define˜ by˜ It follows that˜ :B ρ,κ →B ρ,κ is also an endomorphism. For the constant we will show that: for allα,β ∈B ρ,κ . Fix α, β ∈ B and fixα,β ∈B such that α = i c +α and β = i c +β. Also fix φ ∈ B c (r c ), ξ ∈ B s (r s ). Let us abbreviate x s = x s (t) = x(t, φ, ξ, α) and y s = y s (t) = x(t, φ, ξ, β). Additionally, let us abbreviate x c = α( (t, φ), x s (t)) and y c = β( (t, φ), y s (t)). We wish to bound the following: By Proposition 6.8 we have x s , y s ≤ r s and by (6.11) we have |x c |, |y c | ≤ r c +ρr s , hence: By (6.25) we have |x c − y c | ≤ ρ x s − y s 1 + α −β Ẽ y s 1 , and thereby Plugging this into the integral, and using Proposition 6.8 and Proposition 6.15 to bound x s 1 and x s − y s 1 respectively, we obtain Integrating, we obtain our final bound: After factoring out α −β Ẽ ξ 1 , we have proved inequality (6.26) for our definition of λ.
If λ < 1, as assumed in the theorem, then˜ is a contraction mapping on L ∞ (B c , E) ∩B ρ,κ with norm · Ẽ . So by the Banach fixed point theorem, there exists a unique mapα * ∈B such that˜ [α * ] =α * . If we then define α * = i c +α * it follows that [α * ] = α * . Thereby α * is a chart for the center-stable manifold over B c (r c ).

Computer-assisted proofs of main theorems
In this section, we show computer-assisted proofs of our main results for the initialboundary value problem (2.1). All computations are carried out on Windows 10, Intel(R) Core(TM) i7-6700K CPU @ 4.00GHz, and MATLAB 2019a with INTLAB -INTerval LABoratory [22] version 11 and Chebfun -numerical computing with functions [20] version 5.7.0. All codes used to produce the results in this section are freely available from [25].

Proof of Theorem 1.1
It is obvious that the solution of (1.2) has a symmetry u(te −iθ , x) = u(te iθ , x) for a real t. Then, to prove Theorem 1.1, it is sufficient to show that the imaginary part of u(z, x) is a non-zero function at a certain point z ∈ R satisfying z B < z, where z B denotes the blow-up time of (1.1). We take a path˜ θ which bypasses the blow-up point z B as shown in Fig. 3. Here, z B ≈ 0.0119 holds [4] under the periodic boundary condition with the initial data u 0 (x) = 50(1 − cos(2π x)).  We divide each segment into 64 steps and, by using our rigorous integrator introduced in Sects. 2-5, analytically continue to the z C point in Fig. 3. From O to z A , we set θ = π/3 in (2.1). After that, we change θ = −π/3 from z A to z C . For each time step J i (i = 1, 2, . . . , 128), we get an approximate solution by using Chebfun [20] as ,k T (t) e ikωx , ω = 2π, x ∈ (0, 1), t ∈ J i (7.1) with N = 25 and n = 13. We also set m = 0 in Theorem 3.3, which determines the size of variational problem. The profiles of numerically computed Re(ū) and Im(ū) are plotted in Fig. 4.
After 128 time stepping (at z C ), we show that the imaginary part of the solution u(z C , x) is a non-zero function. To prove this, we use the 1 norm of the Fourier coefficients, say Let the solution u(z C , x) of (1.2) at z C and its numerically computed solution be denoted by respectively. We also denote two bi-infinite complex-valued sequences a z C = (a z C k ) k∈Z andā z C = (. . . , 0,ā z C −N , . . . ,ā z C N , 0 . . . ). Then, the imaginary part of the solution at z C follows Im (u(z C , ·)) = Im (ū(z C , ·)) + Im (u(z C , ·)) − Im (ū(z C , ·)) ≥ Im (ū(z C , ·)) − Im (u(z C , ·) −ū(z C , ·)) where ε z C is the point-wise error at z C point and our rigorous computation yields ε z C = 0.5765. Similary, we have an upper bound of the imaginary part given by Im (u(z C , ·)) ≤ Im (ū(z C , ·)) +ε z C . Furthermore, the imaginary part ofū(z C , x) is presented by This is followed by the following fact: for each |k| ≤ N . Summing up for k, we have We rigorously compute (7.2) based on interval arithmetic and the following is given Im (u(z C , ·)) ∈ [660.4935, 661.6465].
This implies Im (u(z C , ·)) > 0. Consequently, we prove that the imaginary part of the solution u(z C , x) is the nonzero function. Then, there exists a branching singularity inside the contour integral in Fig. 3, i.e., {z ∈ C : Re(z) ≤ z C , 0 ≤ Im(z) ≤ 0.00725 √ 3} with z C = 0.0145. The execute time to prove the branching singularity was almost 105.95 s. The results of analytical continuation is shown in Fig. 5. As shown in this figure, the error bounds  Fig. 4b. Lastly, to obtain a lower bound on the location of the blow-up point z B , it suffices to rigorously integrate our initial data u 0 (x) in purely real time with θ = 0 in (2.1) for as long as possible. Again using Chebfun [20] we obtain an approximate solution as in (7.1) using computational parameters N = 20, n = 15, m = 0, and adaptive time stepping. By using our rigorous integrator introduced in Sects. 2-5, we are able to validate our solution until t = 0.0116, hence the blow-up point must satisfy z B ≥ 0.0116. The execute time to give the above lower bound was almost 69.361 s.

Remark 7.1
Generally speaking, the propagation of error estimates makes rigorous integration difficult. As such, both the length of a contour and how close it approaches the blow-up point will affect the results of our rigorous integrator. We believe the longer contour is a principal reason for explaining why we are able to get 0.0003 close to the blowup point from below, but only 0.0031 close from above. The upper bound could also be improved by adjusting the step size more carefully.

Proof of Theorem 1.2
We show the proof of global existence on the straight path θ . To prove the global existence of the solution of (1.2), we check the hypothesis of Theorem 6.2 by using rigorous numerics. More precisely, at t = t i (after i th time stepping), we rigorously have where a (s) (t) = (a k (t)) k∈Z, k =0 andā (s) (t) = (ā k (t)) |k|≤N , k =0 . To construct a trapping region U as in Corollary 6.4 which might contain a(t i ), we fix mildly inflated radii constants r c def = |ā 0 (t i )| + ε i + 0.02 ā (s) (t i ) 1 + ε i r s def = ā (s) (t i ) 1 + ε i .
Then, for μ = ω 2 cos θ , we compute the nearly optimal bounds of ρ explicitly given in Remark 6.5. If such a (positive) ρ is not obtained, then we are unable validate the center stable manifold in a region large enough to contain a(t i ). In such a case, we continue rigorous integration to the next time step, i.e., J i+1 . Figure 6 shows verified results of our rigorous integration for θ = π/3, π/4, π/6 and π/12. Figure 6a shows results of our rigorous integrator in the case of θ = π/3. After 85 steps of the rigorous integration (at t = 0.2125), the hypothesis of Theorem 6.2 holds for r c = 4.741, r s = 6.545 × 10 −3 , ρ = 1.719 × 10 −2 , 1 = 0.9251, and λ = 0.9251, and a(t) is contained in the region U described in Corollary 6.4. We succeed in proving the global existence of solution on π/3 . The execute time was almost 263.2 sec. Figure 6b also shows results of our rigorous integrator in the case of θ = π/4. At t = 0.1575, the hypothesis of Theorem 6.2 holds for r c = 6.588, r s = 1.033 × 10 −2 , ρ = 1.341 × 10 −2 , 1 = 0.8947, and λ = 0.8947, and a(t) is contained in the region U described in Corollary 6.4. Number of time stepping is 63. We also succeed in proving the global existence of solution on π/4 . It took almost 195.3 s to complete the proof.
In the case of θ = π/6, our rigorous integrator with the same setting above fails to prove local inclusion because the peak of solution becomes large ( ā X ≈ 240). So we change the value of m as m = 0 (0 ≤ t ≤ 0.02) and m = 2 (t ≥ 0.02). Then, as shown in Fig. 6c, our rigorous integrator succeeds in including the solution until t = 0.1325 (53 steps). At that time, the hypothesis of Theorem 6.2 holds for r c = 8.138, r s = 1.382 × 10 −2 , ρ = 1.722 × 10 −2 , 1 = 0.9095, and λ = 0.9095, and a(t) is contained in the region U described in Corollary 6.4. Proof of the global existence on π/6 is complete. The execute time was almost 163.2 sec.
Finally, the case θ = π/12 is slightly difficult to complete. The value of ā X is almost 450 because this path is close to the blow-up point. We adapt both the step size of time stepping and the value of m as Then, after 66 time stepping, the hypothesis of Theorem 6.2 holds at t = 0.125 for r c = 8.939, r s = 2.226 × 10 −2 , ρ = 1.914 × 10 −2 , 1 = 0.8838, and λ = 0.8838, and a(t) is contained in the region U described in Corollary 6.4. On π/12 , we also prove the global existence of solution of (2.1). The results of rigorous integrator is shown in Fig. 6d. It took almost 206.9 s to complete the proof.

Conclusion
In this paper, we introduced a computational method for computing rigorous local inclusions of solutions of the Cauchy problems for the nonlinear heat equation ( 2). Afterwards, in Sect. 6, we introduced an approach based on the Lyapunov-Perron method for calculating part of a center-stable manifold, which is then used to prove that an open set of solutions of the given Cauchy problem converge to zero, hence yielding the global existence in time of the solutions. This method is used to prove our second main result, namely Theorem 1.2. The results are presented in details in Sect. 7. Finally, we conclude this paper by discussing some potential extensions and remained problems for the complex-valued nonlinear heat equation (1.2). We believe that our methodology of rigorous numerics could be extended to more general timedependent PDEs. It is interesting to rigorously integrate the solution orbit of high (space) dimensional PDEs. Furthermore, the case where the time is pure imaginary offers us an interesting problem. The computation in [4] suggests strongly that the solution exists globally in time and decays slowly toward the zero solution. Since the semigroup generator becomes conservative in that case, our method in this paper, which depends on the existence of dissipation, cannot be used.