2D Swift–Hohenberg Equation Dynamics

• The two-dimensional Swift–Hohenberg equation is a fourth-order pattern-forming PDE that generates diverse structures such as stripes, hexagons, holes, and localized states.
• It combines narrow-band linear instability with nonlinear saturation in both conserved and nonconserved frameworks to capture fluid-to-crystal transitions and defect dynamics.
• The model underpins advanced numerical methods that preserve energy decay and mass conservation, enabling accurate simulations of complex pattern coexistence and evolution.

The two-dimensional Swift–Hohenberg equation is a fourth-order pattern-forming partial differential equation for a scalar field on $\mathbb{R}^2$, a periodic rectangle, or related two-dimensional geometries. In its canonical semilinear form it is written as

$u_t = r u - (1+\Delta)^2 u - f'(u),$

with the cubic example

$u_t = r u - (1+2\Delta+\Delta^2)u - u^3,$

while an important conserved variant, often identified with the phase-field crystal model, is

$$\partial_t \phi = \alpha\,\nabla^2\!\left[r\,\phi + (\nabla^2+q^2)^2\phi + \phi^3\right].$$

Across these forms, the equation combines a narrow-band linear instability with nonlinear saturation, and in two dimensions it supports stripes, hexagonal states, holes, localized structures, interfaces, and defect-mediated dynamics. Depending on the variant, it is an $L^2$ gradient flow or a conserved gradient flow, and much of its modern theory is organized around onset, amplitude reductions, coexistence, snaking, secondary instabilities, and structure-preserving numerics (Yonekura et al., 30 Jan 2026, Thiele et al., 2013).

1. Governing formulations and variational structure

The nonconserved two-dimensional Swift–Hohenberg equation is commonly studied on rectangles or periodic tori, with $\Delta=\partial_x^2+\partial_y^2$. One form analyzed in detail is

$$u_t(x,y,t) = -\big[u(x,y,t)^3 + (1-\eta)u(x,y,t) + \varepsilon \Delta^2 u(x,y,t) + 2\Delta u(x,y,t)\big],$$

with either homogeneous Neumann conditions

$$\partial_n u = 0,\qquad \partial_n(\Delta u)=0\quad\text{on }\partial\Omega,$$

or periodic boundary conditions. Its continuous energy is

$$E[u] = \int_\Omega \left(\frac14 u^4 + \frac{1-\eta}{2}u^2 + \frac{\varepsilon}{2}(\Delta u)^2 - |\nabla u|^2\right)\,dx\,dy,$$

and it satisfies the dissipation law

$$\frac{dE}{dt} = -\int_\Omega (u_t)^2\,dx\,dy \le 0.$$

This places the nonconserved equation within the class of fourth-order gradient flows whose dynamics decrease a Lyapunov functional (Yonekura et al., 30 Jan 2026).

The conserved Swift–Hohenberg equation replaces $u_t = r u - (1+\Delta)^2 u - f'(u),$0 relaxation by mass-conserving dynamics. With order parameter $u_t = r u - (1+\Delta)^2 u - f'(u),$1, mobility $u_t = r u - (1+\Delta)^2 u - f'(u),$2, and Helmholtz free energy

$u_t = r u - (1+\Delta)^2 u - f'(u),$3

the evolution is

$u_t = r u - (1+\Delta)^2 u - f'(u),$4

Hence

$u_t = r u - (1+\Delta)^2 u - f'(u),$5

Under periodic or no-flux conditions, the spatial integral of $u_t = r u - (1+\Delta)^2 u - f'(u),$6 is conserved, so the mean order parameter

$u_t = r u - (1+\Delta)^2 u - f'(u),$7

is fixed by the initial condition. For steady states on periodic domains, $u_t = r u - (1+\Delta)^2 u - f'(u),$8 is constant and the steady equation is

$u_t = r u - (1+\Delta)^2 u - f'(u),$9

This conserved formulation is used as the simplest microscopic dynamical model of the fluid-to-crystal transition in the phase-field crystal interpretation (Thiele et al., 2013).

Several two-dimensional generalizations alter symmetry or spectral structure without changing the basic role of fourth-order selection. A generalized Swift–Hohenberg equation with quadratic–cubic nonlinearity,

$u_t = r u - (1+2\Delta+\Delta^2)u - u^3,$0

introduces broken reflection symmetry through the $u_t = r u - (1+2\Delta+\Delta^2)u - u^3,$1 term. An anisotropic stochastic form,

$u_t = r u - (1+2\Delta+\Delta^2)u - u^3,$2

penalizes deviations from criticality differently in the $u_t = r u - (1+2\Delta+\Delta^2)u - u^3,$3 and $u_t = r u - (1+2\Delta+\Delta^2)u - u^3,$4 directions. A dispersive form,

$u_t = r u - (1+2\Delta+\Delta^2)u - u^3,$5

adds a linear dispersive term $u_t = r u - (1+2\Delta+\Delta^2)u - u^3,$6, breaking $u_t = r u - (1+2\Delta+\Delta^2)u - u^3,$7 symmetry and inducing propagation along $u_t = r u - (1+2\Delta+\Delta^2)u - u^3,$8 (Chae et al., 30 Nov 2025, Fukuizumi et al., 2022, Balch et al., 2022).

2. Linear instability, onset, and amplitude reductions

For the standard Turing onset on $u_t = r u - (1+2\Delta+\Delta^2)u - u^3,$9, the linear operator

$$\partial_t \phi = \alpha\,\nabla^2\!\left[r\,\phi + (\nabla^2+q^2)^2\phi + \phi^3\right].$$0

has dispersion relation

$$\partial_t \phi = \alpha\,\nabla^2\!\left[r\,\phi + (\nabla^2+q^2)^2\phi + \phi^3\right].$$1

for Fourier mode $$\partial_t \phi = \alpha\,\nabla^2\!\left[r\,\phi + (\nabla^2+q^2)^2\phi + \phi^3\right].$$2. The threshold is $$\partial_t \phi = \alpha\,\nabla^2\!\left[r\,\phi + (\nabla^2+q^2)^2\phi + \phi^3\right].$$3 with critical wavenumber $$\partial_t \phi = \alpha\,\nabla^2\!\left[r\,\phi + (\nabla^2+q^2)^2\phi + \phi^3\right].$$4, and for $$\partial_t \phi = \alpha\,\nabla^2\!\left[r\,\phi + (\nabla^2+q^2)^2\phi + \phi^3\right].$$5 the unstable band is concentrated near $$\partial_t \phi = \alpha\,\nabla^2\!\left[r\,\phi + (\nabla^2+q^2)^2\phi + \phi^3\right].$$6. In a co-moving frame $$\partial_t \phi = \alpha\,\nabla^2\!\left[r\,\phi + (\nabla^2+q^2)^2\phi + \phi^3\right].$$7, the corresponding dispersion relation becomes

$$\partial_t \phi = \alpha\,\nabla^2\!\left[r\,\phi + (\nabla^2+q^2)^2\phi + \phi^3\right].$$8

which is the starting point for planar interface analysis (Hilder et al., 10 Apr 2026).

For the conserved model, linearization about a homogeneous state $$\partial_t \phi = \alpha\,\nabla^2\!\left[r\,\phi + (\nabla^2+q^2)^2\phi + \phi^3\right].$$9 with perturbation $L^2$0 yields

$L^2$1

The homogeneous state loses stability when

$L^2$2

and the most unstable band is centered near $L^2$3. On finite periodic domains the allowed $L^2$4 are discrete, so onset is determined by the lattice of admissible Fourier modes (Thiele et al., 2013).

Near onset, amplitude reductions organize the two-dimensional patterned states. For hexagonal resonance with

$L^2$5

the formal $L^2$6-equivariant amplitude system is

$L^2$7

Within this system, one-mode roll states and triad hexagons appear as equilibria, and traveling-wave reductions along a direction $L^2$8 replace $L^2$9 by derivatives in $\Delta=\partial_x^2+\partial_y^2$0 (Hilder et al., 10 Apr 2026).

Other singular limits produce different envelope equations. In the anisotropic stochastic equation on a large torus near threshold, the modulation ansatz

$\Delta=\partial_x^2+\partial_y^2$1

leads to the stochastic Ginzburg–Landau equation

$\Delta=\partial_x^2+\partial_y^2$2

with anisotropic diffusion coefficients $\Delta=\partial_x^2+\partial_y^2$3 and $\Delta=\partial_x^2+\partial_y^2$4 in $\Delta=\partial_x^2+\partial_y^2$5 and $\Delta=\partial_x^2+\partial_y^2$6 (Fukuizumi et al., 2022). In the quartic-degeneracy regime, the generalized equation

$\Delta=\partial_x^2+\partial_y^2$7

has neutral curve

$\Delta=\partial_x^2+\partial_y^2$8

and at $\Delta=\partial_x^2+\partial_y^2$9 its minimum at $$u_t(x,y,t) = -\big[u(x,y,t)^3 + (1-\eta)u(x,y,t) + \varepsilon \Delta^2 u(x,y,t) + 2\Delta u(x,y,t)\big],$$0 is quartic,

$$u_t(x,y,t) = -\big[u(x,y,t)^3 + (1-\eta)u(x,y,t) + \varepsilon \Delta^2 u(x,y,t) + 2\Delta u(x,y,t)\big],$$1

The corresponding amplitude equation is

$$u_t(x,y,t) = -\big[u(x,y,t)^3 + (1-\eta)u(x,y,t) + \varepsilon \Delta^2 u(x,y,t) + 2\Delta u(x,y,t)\big],$$2

which reduces, when $$u_t(x,y,t) = -\big[u(x,y,t)^3 + (1-\eta)u(x,y,t) + \varepsilon \Delta^2 u(x,y,t) + 2\Delta u(x,y,t)\big],$$3 and $$u_t(x,y,t) = -\big[u(x,y,t)^3 + (1-\eta)u(x,y,t) + \varepsilon \Delta^2 u(x,y,t) + 2\Delta u(x,y,t)\big],$$4, to a complex Swift–Hohenberg equation with real coefficients (Bentley et al., 2020).

3. Extended patterns, coexistence, and localization

In the two-dimensional conserved Swift–Hohenberg equation with $$u_t(x,y,t) = -\big[u(x,y,t)^3 + (1-\eta)u(x,y,t) + \varepsilon \Delta^2 u(x,y,t) + 2\Delta u(x,y,t)\big],$$5, the phase diagram contains three principal spatially extended patterned phases: stripes, bumps, and holes. Here “stripes” are lamellae or rolls, “bumps” are hexagonally coordinated density maxima, and “holes” are hexagonally coordinated density minima. The same diagram contains several thermodynamic coexistence regions between bumps and the uniform liquid, between bumps and stripes, between stripes and holes, and between holes and the uniform state; the linear stability limit of the uniform state is the spinodal (Thiele et al., 2013).

Direct numerical simulations at $$u_t(x,y,t) = -\big[u(x,y,t)^3 + (1-\eta)u(x,y,t) + \varepsilon \Delta^2 u(x,y,t) + 2\Delta u(x,y,t)\big],$$6, $$u_t(x,y,t) = -\big[u(x,y,t)^3 + (1-\eta)u(x,y,t) + \varepsilon \Delta^2 u(x,y,t) + 2\Delta u(x,y,t)\big],$$7, with periodic boundary conditions on $$u_t(x,y,t) = -\big[u(x,y,t)^3 + (1-\eta)u(x,y,t) + \varepsilon \Delta^2 u(x,y,t) + 2\Delta u(x,y,t)\big],$$8 domains and random initial conditions show a sequence of morphology changes as $$u_t(x,y,t) = -\big[u(x,y,t)^3 + (1-\eta)u(x,y,t) + \varepsilon \Delta^2 u(x,y,t) + 2\Delta u(x,y,t)\big],$$9 varies. For small $$\partial_n u = 0,\qquad \partial_n(\Delta u)=0\quad\text{on }\partial\Omega,$$0, labyrinthine stripe states fill the domain. As $$\partial_n u = 0,\qquad \partial_n(\Delta u)=0\quad\text{on }\partial\Omega,$$1 increases, stripes pinch off locally into bumps, often starting at stripe ends or high-curvature regions. Around $$\partial_n u = 0,\qquad \partial_n(\Delta u)=0\quad\text{on }\partial\Omega,$$2 to $$\partial_n u = 0,\qquad \partial_n(\Delta u)=0\quad\text{on }\partial\Omega,$$3, the bump phase dominates in a hexagonal-crystal-like state with numerous defects. At larger $$\partial_n u = 0,\qquad \partial_n(\Delta u)=0\quad\text{on }\partial\Omega,$$4, vacancies appear within the bump crystal, and eventually the crystal “melts” into individual bumps or smaller clusters. As $$\partial_n u = 0,\qquad \partial_n(\Delta u)=0\quad\text{on }\partial\Omega,$$5 decreases from dense bump arrays to disconnected clusters, a sequence of localized states appears: near $$\partial_n u = 0,\qquad \partial_n(\Delta u)=0\quad\text{on }\partial\Omega,$$6 the bump network percolates, and below this the fragmented “solid-in-liquid” inclusions shrink as the solid fraction falls (Thiele et al., 2013).

These localized states are finite patches of one patterned phase embedded in a background of another, or in the uniform phase: clusters of bumps within stripes or liquid, or holes within a bump crystal. Their thermodynamic organization is controlled by coexistence. In the conserved setting, coexistence between phases occurs on binodal curves where $$\partial_n u = 0,\qquad \partial_n(\Delta u)=0\quad\text{on }\partial\Omega,$$7, $$\partial_n u = 0,\qquad \partial_n(\Delta u)=0\quad\text{on }\partial\Omega,$$8, and the pressure $$\partial_n u = 0,\qquad \partial_n(\Delta u)=0\quad\text{on }\partial\Omega,$$9 match between phases. When localized-state branches are plotted against the conserved mean $$E[u] = \int_\Omega \left(\frac14 u^4 + \frac{1-\eta}{2}u^2 + \frac{\varepsilon}{2}(\Delta u)^2 - |\nabla u|^2\right)\,dx\,dy,$$0 at fixed $$E[u] = \int_\Omega \left(\frac14 u^4 + \frac{1-\eta}{2}u^2 + \frac{\varepsilon}{2}(\Delta u)^2 - |\nabla u|^2\right)\,dx\,dy,$$1, they exhibit slanted homoclinic snaking because the addition of new cells requires global mass redistribution. When the response is instead plotted as $$E[u] = \int_\Omega \left(\frac14 u^4 + \frac{1-\eta}{2}u^2 + \frac{\varepsilon}{2}(\Delta u)^2 - |\nabla u|^2\right)\,dx\,dy,$$2 versus the chemical potential $$E[u] = \int_\Omega \left(\frac14 u^4 + \frac{1-\eta}{2}u^2 + \frac{\varepsilon}{2}(\Delta u)^2 - |\nabla u|^2\right)\,dx\,dy,$$3, the slanted snakes straighten into the vertically aligned homoclinic snaking familiar from nonconserved problems. In two dimensions this same thermodynamic effect appears as plateaus in $$E[u] = \int_\Omega \left(\frac14 u^4 + \frac{1-\eta}{2}u^2 + \frac{\varepsilon}{2}(\Delta u)^2 - |\nabla u|^2\right)\,dx\,dy,$$4 and $$E[u] = \int_\Omega \left(\frac14 u^4 + \frac{1-\eta}{2}u^2 + \frac{\varepsilon}{2}(\Delta u)^2 - |\nabla u|^2\right)\,dx\,dy,$$5 across coexistence intervals; the residual slope is a finite-size effect due to interfacial free-energy contributions scaling as $$E[u] = \int_\Omega \left(\frac14 u^4 + \frac{1-\eta}{2}u^2 + \frac{\varepsilon}{2}(\Delta u)^2 - |\nabla u|^2\right)\,dx\,dy,$$6 rather than bulk $$E[u] = \int_\Omega \left(\frac14 u^4 + \frac{1-\eta}{2}u^2 + \frac{\varepsilon}{2}(\Delta u)^2 - |\nabla u|^2\right)\,dx\,dy,$$7 (Thiele et al., 2013).

A distinct localization mechanism arises in the quartic-minimum generalized equation, where two nearby minima of the neutral curve allow one wavelength to be embedded in a background of another wavelength. In that setting, localized states are “wavelength-in-wavelength” structures rather than patches in a homogeneous background, and the amplitude equation supports fronts and localized combinations of two detunings. The same analysis reports that, without the Proctor term, the branches are disconnected, whereas with a nonzero Proctor term they join and display more classical snaking (Bentley et al., 2020).

4. Secondary instabilities, defect dynamics, and defect-mediated transport

Two-dimensional stripe states are not generically stable. In the generalized equation

$$E[u] = \int_\Omega \left(\frac14 u^4 + \frac{1-\eta}{2}u^2 + \frac{\varepsilon}{2}(\Delta u)^2 - |\nabla u|^2\right)\,dx\,dy,$$8

with $$E[u] = \int_\Omega \left(\frac14 u^4 + \frac{1-\eta}{2}u^2 + \frac{\varepsilon}{2}(\Delta u)^2 - |\nabla u|^2\right)\,dx\,dy,$$9 and $$\frac{dE}{dt} = -\int_\Omega (u_t)^2\,dx\,dy \le 0.$$0, rolls are steady $$\frac{dE}{dt} = -\int_\Omega (u_t)^2\,dx\,dy \le 0.$$1-periodic solutions depending only on $$\frac{dE}{dt} = -\int_\Omega (u_t)^2\,dx\,dy \le 0.$$2. Linearization about a roll leads, after Bloch–Floquet reduction, to the family of self-adjoint operators

$$\frac{dE}{dt} = -\int_\Omega (u_t)^2\,dx\,dy \le 0.$$3

For sufficiently small $$\frac{dE}{dt} = -\int_\Omega (u_t)^2\,dx\,dy \le 0.$$4 and all $$\frac{dE}{dt} = -\int_\Omega (u_t)^2\,dx\,dy \le 0.$$5, the roll is spectrally unstable; unstable Bloch parameters lie in a band of width $$\frac{dE}{dt} = -\int_\Omega (u_t)^2\,dx\,dy \le 0.$$6 around the resonant set, and the paper proves a transition from spectral to nonlinear instability by constructing perturbations localized near the maximally unstable Bloch mode. The resulting nonlinear theorem yields an exit time

$$\frac{dE}{dt} = -\int_\Omega (u_t)^2\,dx\,dy \le 0.$$7

for perturbations of size $$\frac{dE}{dt} = -\int_\Omega (u_t)^2\,dx\,dy \le 0.$$8 in $$\frac{dE}{dt} = -\int_\Omega (u_t)^2\,dx\,dy \le 0.$$9, with deviation at time $u_t = r u - (1+\Delta)^2 u - f'(u),$00 measured in $u_t = r u - (1+\Delta)^2 u - f'(u),$01 (Chae et al., 30 Nov 2025).

A different instability theory arises when stripe patterns are coupled to a mean flow. In two generalized Swift–Hohenberg models with a projection operator chosen so that single-mode stripes are exact solutions, stripes become unstable not only to the usual Eckhaus and zigzag modes but also to skew-varicose, oscillatory skew-varicose, and cross-roll instabilities. The exact stripe solution has the form

$u_t = r u - (1+\Delta)^2 u - f'(u),$02

For stress-free boundary conditions, the skew-varicose boundary is

$u_t = r u - (1+\Delta)^2 u - f'(u),$03

and the oscillatory skew-varicose boundary, present in one of the mean-flow models, is

$u_t = r u - (1+\Delta)^2 u - f'(u),$04

for small $u_t = r u - (1+\Delta)^2 u - f'(u),$05. At large coupling, the cross-roll instability can eliminate the stable stripe region completely (Weliwita et al., 2011).

The addition of linear dispersion produces yet another two-dimensional defect class. In the dispersive Swift–Hohenberg equation

$u_t = r u - (1+\Delta)^2 u - f'(u),$06

the term $u_t = r u - (1+\Delta)^2 u - f'(u),$07 generates traveling ripples and spatially extended dislocations called seam defects. Near threshold, the amplitude equation is

$u_t = r u - (1+\Delta)^2 u - f'(u),$08

a special case of the anisotropic complex Ginzburg–Landau equation. In this reduction, seams correspond to spiral waves of the envelope, and chains of spiral cores satisfy the velocity and spacing laws

$u_t = r u - (1+\Delta)^2 u - f'(u),$09

In the strong-dispersion limit for the one-dimensional traveling-wave reduction, the amplitude–wavenumber and velocity relations are

$u_t = r u - (1+\Delta)^2 u - f'(u),$10

with existence restricted to $u_t = r u - (1+\Delta)^2 u - f'(u),$11 (Balch et al., 2022).

5. Rigorous bifurcation, interfaces, heterogeneity, and symmetry

Close to a Turing instability, two-dimensional Swift–Hohenberg-type equations admit slow planar invasion fronts in general directions. For

$u_t = r u - (1+\Delta)^2 u - f'(u),$12

with $u_t = r u - (1+\Delta)^2 u - f'(u),$13 symmetry and a smallness condition on quadratic resonant triads, a spatial-dynamics and non-standard center-manifold analysis gives a finite-dimensional locally invariant center manifold for angles satisfying

$u_t = r u - (1+\Delta)^2 u - f'(u),$14

On this manifold, the reduced traveling-wave equations are the anisotropic hexagonal Ginzburg–Landau system

$u_t = r u - (1+\Delta)^2 u - f'(u),$15

and, for $u_t = r u - (1+\Delta)^2 u - f'(u),$16, $u_t = r u - (1+\Delta)^2 u - f'(u),$17, $u_t = r u - (1+\Delta)^2 u - f'(u),$18, and $u_t = r u - (1+\Delta)^2 u - f'(u),$19, persistent heteroclinic connections from stripes or down-hexagons to the trivial state lift to slow-moving planar pattern interfaces of the original PDE,

$u_t = r u - (1+\Delta)^2 u - f'(u),$20

The front speed scales as $u_t = r u - (1+\Delta)^2 u - f'(u),$21, and the anisotropic coefficients $u_t = r u - (1+\Delta)^2 u - f'(u),$22 make the interface strongly direction-dependent (Hilder et al., 10 Apr 2026).

Compact heterogeneity provides a separate rigorous route to fully localized two-dimensional states. For

$u_t = r u - (1+\Delta)^2 u - f'(u),$23

the trivial state is linearly unstable only inside the disc $u_t = r u - (1+\Delta)^2 u - f'(u),$24 and linearly stable outside. The essential spectrum of the linearization about a steady state is

$u_t = r u - (1+\Delta)^2 u - f'(u),$25

and separation in polar coordinates yields localized eigenfunctions

$u_t = r u - (1+\Delta)^2 u - f'(u),$26

Under a one-dimensional kernel hypothesis, analytic Crandall–Rabinowitz bifurcation gives local branches in the symmetry-restricted spaces $u_t = r u - (1+\Delta)^2 u - f'(u),$27, with $u_t = r u - (1+\Delta)^2 u - f'(u),$28 producing axisymmetric spots and $u_t = r u - (1+\Delta)^2 u - f'(u),$29 dipole-type states. For $u_t = r u - (1+\Delta)^2 u - f'(u),$30 the local bifurcation is a pitchfork; for $u_t = r u - (1+\Delta)^2 u - f'(u),$31 it is transcritical when

$u_t = r u - (1+\Delta)^2 u - f'(u),$32

and a pitchfork otherwise. Analytic global continuation then extends these branches to large amplitude. The primary bifurcating branch alternates between an axisymmetric spot and a non-axisymmetric dipole, depending on the width $u_t = r u - (1+\Delta)^2 u - f'(u),$33 of the heterogeneity (Hill et al., 3 Oct 2025).

Computer-assisted analysis has also validated periodic two-dimensional branches with prescribed space-group symmetry. For the steady equation

$u_t = r u - (1+\Delta)^2 u - f'(u),$34

on the hexagonal lattice, a Fourier representation on the transformed reciprocal lattice combined with Newton–Kantorovich radii-polynomial estimates proves the existence of smooth periodic solutions with $u_t = r u - (1+\Delta)^2 u - f'(u),$35 and $u_t = r u - (1+\Delta)^2 u - f'(u),$36 symmetry, as well as entire branches parameterized by Chebyshev series. The global tiling is generated by two symmetric triangles in the $u_t = r u - (1+\Delta)^2 u - f'(u),$37 case and by a single symmetric hexagon in the $u_t = r u - (1+\Delta)^2 u - f'(u),$38 case (Blanco, 13 Feb 2026).

A geometrically distinct extension replaces the Euclidean plane by the Poincaré disc. There the equation

$u_t = r u - (1+\Delta)^2 u - f'(u),$39

supports $u_t = r u - (1+\Delta)^2 u - f'(u),$40-periodic H-planforms on hyperbolic lattices, radial localized states, and horocyclic traveling waves, with the hyperbolic spectral shift

$u_t = r u - (1+\Delta)^2 u - f'(u),$41

replacing the Euclidean $u_t = r u - (1+\Delta)^2 u - f'(u),$42 dispersion (Chossat et al., 2013).

6. Structure-preserving numerics and validated computation

Because the two-dimensional Swift–Hohenberg equation is a fourth-order gradient flow, numerical schemes are often designed to preserve a discrete energy law, guarantee unique solvability, and remain stable for large time steps. Several families of methods do this while targeting different variants and boundary conditions.

Scheme	Core formulation	Established properties
Unconditionally energy-stable isogeometric DNS for conserved SH/PFC	Isogeometric Analysis with NURBS in space and a second-order perturbation of the trapezoidal rule in time	Free-energy-decreasing discrete solutions for any time step and mesh size; used for periodic $u_t = r u - (1+\Delta)^2 u - f'(u),$43 and $u_t = r u - (1+\Delta)^2 u - f'(u),$44 domains (Thiele et al., 2013)
Smooth–convex–concave splitting finite differences	Linearly implicit treatment of the 2D equation $u_t = r u - (1+\Delta)^2 u - f'(u),$45	Unique solvability; discrete energy decay under the stated $u_t = r u - (1+\Delta)^2 u - f'(u),$46 condition or unconditionally if $u_t = r u - (1+\Delta)^2 u - f'(u),$47; boundedness of numerical solutions; first-order time and second-order space error estimate (Yonekura et al., 30 Jan 2026)
Stabilizing Correction finite differences for SH3 and SH35	Semi-implicit midpoint discretization with operator splitting and positive-definite prefactor $u_t = r u - (1+\Delta)^2 u - f'(u),$48	Unconditional stability verified numerically; strict implementation of the Lyapunov functional; second-order accuracy in time and space; supports GDBC and PBC, ramps and Gaussian forcing (Coelho et al., 2020)
Fully discrete BDF3 finite differences	Fully implicit BDF3 in time and second-order finite differences in space for $u_t = r u - (1+\Delta)^2 u - f'(u),$49	Uniquely solvable; unconditionally energy stable in a modified energy; bounded in maximum norm; $u_t = r u - (1+\Delta)^2 u - f'(u),$50-error estimate of order $u_t = r u - (1+\Delta)^2 u - f'(u),$51 under the stated time-step condition (Zhao et al., 2023)
IEQ–DG methods	Mixed DG spatial discretization plus Invariant Energy Quadratization time stepping	Linear fully discrete schemes; unconditionally energy stable; free-energy-decaying discrete solutions irrespective of time step and mesh size; comparable high-accuracy 2D patterns (Liu et al., 2019)

In the conserved setting, direct numerical simulation has been used not only as a time integrator but as a phase-diagram tool. The two-dimensional cSH/PFC computations employed periodic domains of size $u_t = r u - (1+\Delta)^2 u - f'(u),$52, random initial conditions, and diagnostics including the mean order parameter

$u_t = r u - (1+\Delta)^2 u - f'(u),$53

the amplitude measure

$u_t = r u - (1+\Delta)^2 u - f'(u),$54

the mean free-energy density relative to the homogeneous state, and the mean grand potential density

$u_t = r u - (1+\Delta)^2 u - f'(u),$55

These diagnostics expose coexistence plateaus in $u_t = r u - (1+\Delta)^2 u - f'(u),$56 and $u_t = r u - (1+\Delta)^2 u - f'(u),$57 and allow localized states in two dimensions to be organized thermodynamically (Thiele et al., 2013).

Nonuniform forcing and boundary effects are numerically important even for the nonconserved equation. In the cubic and quintic two-dimensional equations with spatially varying control parameter,

$u_t = r u - (1+\Delta)^2 u - f'(u),$58

the discrete Lyapunov functional decreases monotonically for all reported test cases. Generalized Dirichlet boundary conditions favor perpendicular incidence of rolls at sidewalls and can increase defect density, whereas periodic boundary conditions reduce boundary-induced defects. Ramps and Gaussian forcings produce localization, target patterns, and wavelength filtering, and in the SH35 regime subcritical forcing supports localized stripes (Coelho et al., 2020).

Across these methods, a persistent theme is that faithful computation of two-dimensional Swift–Hohenberg dynamics is tied to the variational or conserved structure of the PDE. The numerical literature therefore emphasizes free-energy decay, exact or approximate mass conservation where appropriate, and discretizations adapted to the fourth-order operator, rather than purely explicit time stepping. This suggests a broad methodological consensus: in two dimensions, reliable computation of coexistence, localization, and defect evolution depends strongly on structure-preserving discretization.