1 Introduction
We consider the AllenCahn equation
(1.1) 
where the unknown . The parameter is called the mobility coefficient and we fix it as a constant for simplicity. The nonlinear term takes the form , where is the standard double well. To minimize technicality, we take the spatial domain in (1.1) as the periodic torus . With some additional work our analysis can be extended to physical dimensions . The system (1.1) arises as a gradient flow of a GinzburgLandau type energy functional , where
(1.2) 
The basic energy conservation law takes the form
(1.3) 
It follows that
(1.4) 
Besides the type conservation law, there is also type control. Due to special form of the nonlinearity, for smooth solutions we have the maximum principle
(1.5) 
As a consequence the longtime wellposedness and regularity is not an issue for (1.1).
The objective of this work is to establish strong stability of a fourth order in time operating splitting algorithm applied to the AllenCahn equation (1.1). This is a part of our ongoing program to develop a new theory for the rigorous analysis of stability and convergence of operatorsplitting methods applied to dissipativetype problems. Due to various subtle technical obstructions, there were very few rigorous results on the analysis of the operatorsplitting type algorithms for the AllenCahn equation, the CahnHilliard equation and similar models. Prior to our recent series of works [17, 18], most existing results in the literature are conditional one way or another. To put things into perspective, we briefly review a few closely related representative works and more recent developments.

The work of GideyReddy. Gidey and Reddy considered in [24] a convective CahnHilliard model of the form
(1.6) where . They adopted an operatorsplitting of (1.6) into the hyperbolic part, nonlinear diffusion part and diffusion part respectively. Several conditional results concerning certain weak solutions were obtained.

The work of WengZhaiFeng. In [25], Weng, Zhai and Feng studied a viscous CahnHilliard model:
(1.7) where the parameter . The authors employed a fast explicit Strangtype operator splitting and proved the stability and the convergence (see Theorem 1 on pp. 7 of [25]) under the assumption that , stay bounded, and satisfy a technical condition . Here denotes the numerical solution.

The work of ChengKurganovQuTang. In [23], Cheng, Kurganov, Qu and Tang considered the CahnHilliard equation
(1.8) and the MBE equation
(1.9) Concerning the CahnHilliard equation, the authors considered a Strangtype splitting approximation of the form
(1.10) where
(1.11) and is the nonlinear propagator
(1.12) Various conditional results were given in [23] but the rigorous analysis of energy stability was a longstanding open problem. This problem and several related open problems were settled in our recent proof [18].
In recent [17], we carried out the first energystability analysis of a first order operatorsplitting approximation of the CahnHilliard equation. More precisely denote as the exact PDE solution to the CahnHilliard equation . We considered a splitting approximation of the form:
(1.13) 
where and solves
(1.14) 
We introduced a novel modified energy and rigorously proved monotonic decay of the new modified energy which is coercive in sense. Moreover we obtained uniform control of higher Sobolev regularity and rigorously justified the first order convergence of the method on any finite time interval.
In [18], we settled the difficult open problem of energy stability of the Strangtype algorithm applied to the CahnHilliard equation which was introduced in the work of Cheng, Kurganov, Qu and Tang [23]. One should note that the new theoretical framework developed in [18] is completely different from the first order case [17]. In the secondorder case, for most generic data one no longer has strict energymonotonicity at disposal and several new ideas such as dichotomy analysis, an absorbingset approach were introduced in [18] in order to settle long time unconditional (i.e. independent of timestep or time interval) Sobolev bounds of the numerical solution.
In this work we develop further the program initiated in [17, 18] and turn to the analysis of higher order in time splitting methods. A folklore issue is that operatorsplitting methods with order higher than two require negative timestepping, i.e. backward time integration for each time step. In [33], Sheng considered dimensional splitting for the twodimensional parabolic problem . Using semidiscretization one obtains the ODE system , where the matrices and correspond to the discretization of and respectively. One is then naturally led to the approximation of . In [33], Sheng showed that if the noncommuting matrices and have eigenvalues in the left half of the complexplane, and one employs the approximation
(1.15)  
(1.16) 
then the highest order of a stable approximation is two even if is chosen to be large. Put it differently, Sheng’s fundamental result states that for an order () partitioned splitstep schemes, at least one of the solution operators must be applied with negative time step, i.e. backwardly. In [34] Suzuki adopted an elegant timeordering principle from quantum mechanics and proved a general nonexistence theorem of positive decomposition for high order splitting approximations. In [35], Goldman and Kaper strengthened these results further and showed that even with partitioned schemes, each solution operator within a convex partition must be performed with at least one negative/backward fractional time step.
For deterministic Hamiltoniantype systems, backward time stepping in general does not create instabilities and highorder operating splitting methods have shown promising effectiveness [32]. On the other hand, due to the negative time stepping, there are some arguments that higherthanthree operator splitting methods cannot be applied to parabolic equations with diffusive terms ([33, 36, 37]). As we shall explain momentarily, these concerns are not unsubstantiated and even turn up in the welldefinedness of these algorithms.
Although rigorous analysis of these issues were not available before, recently Cervi and Spiteri [31] considered three thirdorder operatingsplitting methods and demonstrated via extensive numerical simulations the effectiveness of higher order splitting methods for representative cardiac electrophysiology simulations. These compelling numerical evidences propel us to reexamine in detail the stability property of general negative/backward timestepping methods in parabolic problems. Indeed, the very purpose of this work is to break these aforementioned technical barriers and establish a new stability theory for these problems.
To set the stage and minimize technicality, we consider the AllenCahn equation (1.1) posed on the periodic torus . To build some intuition, let and consider
(1.17) 
and let solve the equation
(1.18) 
An explicit formula for is readily available, indeed it is not difficult to check that
(1.19) 
Define a secondorder Strangtype propagator
(1.20) 
Following Yoshida [29], we consider a 4th order integrator obtained by a symmetric repetition (product) of the 2nd order integrator:
(1.21) 
where
(1.22) 
Somewhat surprisingly, we first show that the propagator is illdefined if one does not introduce a judiciously chosen spectral cutoff.
Proposition 1.1 (Illdefinedness of with a spectral cutoff).
The propagator is illdefined in general.
The proof of Proposition 1.1 is given in Section 2. Armed with this important observation, we are led to introduce a spectral cutoff condition for the propagators. Let be an integer and we shall consider the projection operator defined for via the relation
(1.23) 
where denotes the Fourier coefficient of . We introduce the following spectral cutoff condition.
Definition (Spectral condition). Let be the spectral truncation parameter. We shall say satisfy the spectral condition if for some .
Remark 1.1.
This constraint in is reminiscent of the CFL condition in hyperbolic problems. Here in the parabolic setting we require that the operator to remain bounded when restricted to the spectral cutoff .
Let . We now consider the following modified propagators:
(1.24) 
Theorem 1.1 (Stability of negativetime splitting methods).
Let , and consider the AllenCahn equation (1.1) on the onedimensional periodic torus . Assume the spectral cutoff condition for some . Assume the initial data ( is an integer). Let and define
(1.25) 
There exists a constant depending only on , and , such that if , then
(1.26) 
where depends on (, , , ).
Remark 1.2.
One can also show the fourthorder convergence of the operator splitting approximation. Namely if and let be the exact PDE solution corresponding to initial data . Let bet the same as in Theorem 1.1. Then for any , we have
(1.27) 
where depends on (, , , ). The regularity assumption on initial data can be lowered. We shall refrain from stating such a pedestrian result here and leave its justification to interested readers as exercises.
The rest of this paper is organized as follows. In Section we set up the notation and collect some preliminary lemmas. In Section we give the proof of Theorem 1.1.
2 Notation and preliminaries
For any two positive quantities and , we shall write or if for some constant whose precise value is unimportant. We shall write if both and hold. We write if the constant depends on some parameter . We shall write if and if .
We shall denote if for some sufficiently small constant . The smallness of the constant is usually clear from the context. The notation is similarly defined. Note that our use of and here is different
from the usual Vinogradov notation in number theory or asymptotic analysis.
For any , we denote , and . Also occasionally we use the Japanese bracket notation:
We denote by the usual periodic torus. For and any function , we denote the Lebesgue norm of as
If is a sequence of complex numbers and is the index set, we denote the discrete norm as
(2.1) 
For example, . If
is a vectorvalued function, we denote
, and . We use similar convention for the corresponding discrete norms for the vectorvalued case.We use the following convention for the Fourier transform pair:
(2.2) 
and denote for ,
(2.3a)  
(2.3b) 
Proof of Proposition 1.1.
To prove Proposition 1.1, it suffices for us to examine the following statement. Consider where is the periodic Dirac comb on . Let and consider
(2.4) 
Claim: for any .
Proof of Claim. Observe that the Fourier coefficients of are all positive and
(2.5) 
Clearly the desired conclusion follows. ∎
3 Proof of Theorem 1.1
Lemma 3.1 (Onestep stability).
Let , and assume the spectral condition for some constant . Suppose , and . Let . There exists sufficiently small such that if , then
(3.1) 
where depends on (, , , , , ).
Remark 3.1.
Define
(3.2)  
(3.3)  
(3.4)  
(3.5) 
Later we shall apply Lemma 3.1 times. In particular we need to bound the expression
(3.6) 
Since the constant depends on the norm of the iterates, it is of importance to give uniform control of the norm. To resolve this issue, we first choose , and such that
(3.7) 
We choose such that
(3.8) 
Clearly in the course of iteration, the norm of the iterates never exceeds . In particular
(3.9) 
Proof of Lemma 3.1.
We begin by noting that for the ODE
(3.10) 
we have the estimate (note that
is an algebra for )(3.11) 
It follows that for sufficiently small,
(3.12)  
(3.13) 
The desired estimates then easily follows from the above using the spectral condition. ∎
Lemma 3.2 (Multistep stability and higher regularity).
Let , and assume the spectral condition for some constant . Suppose and for some constant . Define
(3.14) 
where , . For , define
(3.15) 
where . There exist and such that if (we may assume ), then
(3.16) 
where depends on (, , ).
Proof of Lemma 3.2.
Observe that
(3.17) 
where . Clearly the operators , , fulfill the conditions of Lemma 3.1. The estimates then follow easily from the computations outlined in Remark 3.1 with some necessary adjustment of the constants.
To establish (3.16), we note that for ,
(3.18) 
provided . We then rewrite as
(3.19) 
where . One can then bootstrap the higher regularity from this. We omit further details. ∎
Remark 3.2.
It is not difficult to check that under the assumption of uniform high Sobolev regularity, we have (below we assume )
(3.20) 
It follows that
(3.21) 
Lemma 3.3 (Control of the energy flux).
Let and . Suppose satisfies and
(3.22) 
Then
(3.23) 
where depends only on . Furthermore if and where is a sufficiently small constant depending on (, ), then
(3.24) 
where depends only on (, ).
Proof.
Lemma 3.4 (Onestep strict energy dissipation with nontrivial energy flux).
Let , and . Suppose with and satisfies
(3.25) 
where is a given constant. There exists sufficiently small such that if , then
(3.26) 
Proof of Lemma 3.4.
By using Remark 3.2, it suffices for us to show
(3.27) 
for some . Denote and observe that
(3.28) 
It follows that for sufficiently small,
(3.29) 
where is a constant. Now we clearly have
(3.30) 
The desired result clearly follows. ∎
Theorem 3.1.
Let , and . Assume and . Define and
(3.31) 
There exists sufficiently small such that if , then
(3.32) 
where depends only on (, , ).
Proof of Theorem 3.1.
Denote
(3.33) 
where the constants , are the same as in Lemma 3.3. In the argument below we shall assume is sufficiently small. The needed smallness (i.e. the existence of ) can be easily worked out by fulfilling the conditions needed in Lemma 3.2, Lemma 3.3 and Lemma 3.4.
By Lemma 3.2, we can find such that
(3.34) 
Claim: We have
(3.35) 
To prove the claim we argue by contradiction. Suppose is the first integer such that
(3.36) 
By Lemma 3.3 we must have
(3.37) 
Since , we have for some integer . By using smoothing estimates we obtain
(3.38) 
where depends on (, , ). Since depends on (, , ), we have depends only on (, , ). By (3.37), (3.38) and Lemma 3.4, we obtain for sufficiently small that
(3.39) 
which is clearly a contradiction to (3.36). Thus we have proved the claim. ∎
References
 [1] A. Baskaran, J. S. Lowengrub, C. Wang, and S. M. Wise. Convergence Analysis of a Second Order Convex Splitting Scheme for the Modified Phase Field Crystal Equation. SIAM J. Numer. Anal., 51(2013), 2851šC2873.
 [2] J.W. Cahn, J.E. Hilliard. Free energy of a nonuniform system. I. Interfacial energy free energy, J. Chem. Phys. 28 (1958) 258–267.
 [3] L.Q. Chen, J. Shen. Applications of semiimplicit Fourierspectral method to phase field equations. Comput. Phys. Comm., 108 (1998), pp. 147–158.
 [4] A. Christlieb, J. Jones, K. Promislow, B. Wetton, M. Willoughby. High accuracy solutions to energy gradient flows from material science models. J. Comput. Phys. 257 (2014), part A, 193–215.
 [5] C.M. Elliott and A.M. Stuart. The global dynamics of discrete semilinear parabolic equations. SIAM J. Numer. Anal., 30 (1993), pp. 1622–1663.
 [6] D. J. Eyre. Unconditionally gradient stable time marching the CahnHilliard equation. Computational and mathematical models of microstructural evolution (San Francisco, CA, 1998), 39–46, Mater. Res. Soc. Sympos. Proc., 529, MRS, Warrendale, PA, 1998.
 [7] X. Feng, T. Tang and J. Yang. Long time numerical simulations for phasefield problems using padaptive spectral deferred correction methods. SIAM J. Sci. Comput. 37 (2015), no. 1, A271–A294.
 [8] H. Gomez and T.J.R. Hughes. Provably unconditionally stable, secondorder timeaccurate, mixed variational methods for phasefield models. J. Comput. Phys., 230 (2011), pp. 53105327
 [9] Z. Guan, C. Wang and S. M. Wise, A convergent convex splitting scheme for the periodic nonlocal CahnHilliard equation, Numer. Math., 128 (2014), 377–406.
 [10] J. Guo, C. Wang, S. M. Wise and X. Yue, An convergence of a secondorder convexsplitting, finite difference scheme for the threedimensional CahnHilliard equation, Commu. Math. Sci., 14 (2016), 489–515
 [11] Y. He, Y. Liu and T. Tang. On large timestepping methods for the CahnHilliard equation. Appl. Numer. Math., 57 (2007), 616–628.
 [12] F. Liu and J. Shen. Stabilized semiimplicit spectral deferred correction methods for AllenCahn and CahnHilliard equations. Math. Methods Appl. Sci. 38 (2015), no. 18, 4564–4575.
 [13] D. Li, Z. Qiao and T. Tang. Characterizing the stabilization size for semiimplicit Fourierspectral method to phase field equations, SIAM J. Numer. Anal., 54 (2016), 1653–1681
 [14] D. Li, F. Wang and K. Yang. An improved gradient bound for 2D MBE. Journal of Differential Equations, 269(12): 1116511171, 2020.
 [15] D. Li, T. Tang. Stability of the SemiImplicit Method for the CahnHilliard Equation with Logarithmic Potentials. Ann. Appl. Math., 37 (2021), 31–60.
 [16] D. Li, C. Quan, and T. Tang. Stability analysis for the implicitexplicit discretization of the CahnHilliard equation.
 [17] D. Li and C. Quan. The operatorsplitting method for CahnHilliard is stable. arXiv:2107.01418, 2021.
 [18] D. Li and C. Quan. On the energy stability of Strangsplitting for CahnHilliard. arXiv:2107.05349, 2021.
 [19] C.B. Schönlieb and A. Bertozzi. Unconditionally stable schemes for higher order inpainting. Commun. Math. Sci. 9 (2011), no. 2, 413–457.
 [20] J. Shen and X. Yang. Numerical approximations of AllenCahn and CahnHilliard equations. Discrete Contin. Dyn. Syst. A, 28 (2010), 1669–1691.
 [21] J. Shen, J. Xu and J. Yang. A new class of efficient and robust energy stable schemes for gradient flows. SIAM Rev. 61 (2019), no. 3, 474–506.
 [22] H. Song and C. Shu. Unconditional energy stability analysis of a second order implicitexplicit local discontinuous Galerkin method for the CahnHilliard equation. J. Sci. Comput. 73 (2017), no. 23, 1178–1203.
 [23] Y. Cheng, A. Kurganov, Z. Qu, and T. Tang. Fast and stable explicit operator splitting methods for phasefield models. Journal of Computational Physics, 303:45–65, 2015.
 [24] H.H. Gidey and B.D. Reddy. Operatorsplitting methods for the 2D convective CahnHilliard equation. Computers & Mathematics with Applications, 2019, 77(12): 3128–3153.

[25]
Z. Weng, S. Zhai and X. Feng. Analysis of the operator splitting scheme for the CahnHilliard equation with a viscosity term. Numer. Meth. for Partial Differential Equations. 35 (2019), no. 6, 1949–1970.
 [26] C. Xu and T. Tang. Stability analysis of large timestepping methods for epitaxial growth models. SIAM J. Numer. Anal. 44 (2006), no. 4, 1759–1779.
 [27] J. Zhu, L.Q. Chen, J. Shen, and V. Tikare. Coarsening kinetics from a variablemobility CahnHilliard equation: Application of a semiimplicit Fourier spectral method, Phys. Rev. E (3), 60 (1999), pp. 3564–3572.
 [28] R. I. Issa. Solution of the implicitly discretised fluid flow equations by operatorsplitting. Journal of Computational Physics 62, no. 1 (1986): 40–65.
 [29] H. Yoshida. Construction of higher order symplectic integrators. Physics letters A 150, no. 5–7 (1990): 262268.
 [30] R. I. McLachlan and G. R. W. Quispel. Splitting methods. Acta Numerica 11 (2002): 341–434.
 [31] J. Cervi and R. J. Spiteri. Highorder operator splitting for the bidomain and monodomain models. SIAM Journal on Scientific Computing 40, no. 2 (2018): A769–A786.
 [32] M. Thalhammer, M. Caliari and C. Neuhauser. Highorder timesplitting Hermite and Fourier spectral methods. J. Comput. Phys. 228 (2009), pp. 822–832.
 [33] Q. Sheng. Solving linear partial differential equations by exponential splitting. IMA J. Numer. Anal. 9 (1989), pp. 199–212.
 [34] M. Suzuki. General theory of fractal path integrals with applications to manybody theories and statistical physics. J. Math. Phys., 32 (1991), pp. 400–407.
 [35] G. Goldman and T.J. Kaper. Nthorder operator splitting schemes and nonreversible systems. SIAM J. Numer. Anal. 33 (1996), pp. 349–367.
 [36] S.A. Chin. Quantum statistical calculations and symplectic corrector algorithms. Phys. Rev. E. 69 (2004), 046118.
Comments
There are no comments yet.