NONCLASSICAL SYMMETRY SOLUTIONS FOR 4 TH ORDER PHASE FIELD REACTION-DIFFUSION

Using a nonclassical symmetry of nonlinear reaction-diffusion equations, 1 some exact multi-dimensional time-dependent solutions are constructed for a fourth-order 2 Allen-Cahn-Hilliard equation. This models a phase field that gives a phenomenological 3 description of a two-phase system near the critical temperature. Solutions are given for 4 the changing phase of a cylindrical or spherical inclusion, allowing for a ’mushy zone’ 5 with mixed state that is controlled by imposing a pure state at the boundary. The diffusion 6 coefficients for transport of one phase through the mixture, depend on the phase field value, 7 since the physical structure of the mixture depends on the relative proportions of the two 8 phases. A source term promotes stability of both of the pure phases but this tendency may 9 be controlled or even reversed through the boundary conditions. 10


Introduction
Material transport problems with phase change are conveniently modelled by Stefan free boundary conditions, in which there is assumed to be a moving sharp interface between phases.In practice, the interface of solidification may consist of a "mushy" zone, within which there is dendritic infiltration of one phase by another.Similarly, crystal dislocations may migrate between one solid crystal phase and another through a region of crystal overlap that may be significant at the microscale (e.g., [1]).Observable transition zones with two-phase mixtures may occur in environments wherein the temperature varies little from its critical value.For multi-phase dynamical modelling, the next level of realism beyond sharp interfaces is the phase field theory, in which a continuously varying scalar field θ(x, t) interpolates Phase 1 (θ = 1) and Phase 0 (θ = 0).θ is a monotonic function of the concentration of Phase 1.It may be regarded as a surrogate for the relative concentration of one phase.As in most continuum representations of mixtures, after averaging over a representative elementary volume that is large compared to individual particles such as microscopic dendrites, quantities such as density and phase concentration that are discontinuous at the microscale are represented as smoothly varying functions.The concept of the phase variable is a mathematical generalisation of the older concept of a physically measurable order parameter [2].The classic example is the magnetisation of a ferromagnetic material.Only since the 1940s have we been able to calculate the critical temperature for the simplest phase transition, which is the magnetisation of the spatially uniform two-dimensional Ising ferromagnet.In practice, phase transitions are not instantaneous and during the transition, the phase structure of a material will be heterogeneous at the microscale.The scalar field theory of a single phase variable, while being a relatively simple model, will still predict some spatially non-uniform phase structures through a transition zone.The transition zone is often referred to as the "mushy" zone, following the analogy of a two-phase water-ice mixture.The morphology of the transition zone is a major determinant of the physical properties of alloys (e.g., [3]).The Cahn and Hilliard phase field model of phase separation [4] showed why the transition zone is stable only at length scales below a small characteristic value.Following the works by Allen and Cahn [5], and by Cahn and Hilliard [4], phase field theories have been studied for several decades [6][7][8][9][10].
The core of usual simplified versions of Allen-Cahn (or of Ginzburg-Landau [11]) is an equation for direct relaxation to the minima of a double-well potential.The resulting partial differential equation (PDE) is a standard reaction-diffusion equation of the Fitzhugh-Nagumo type, where β > 0 and s are constants.Here, the two stable steady states are the pure phases θ = 0 and θ = 1.Without a loss of generality, s > 0; otherwise we could interchange the labels of the two phases.
The cubic source term drives the dynamics towards bi-stable fixed points, with the zones of attraction separated by the unstable fixed-point at θ = 0.5.The diffusion term allows for forward transport of existing Phase 1 material through the mixture, opposite to the direction of its own concentration gradient.Equivalently, Phase 0 diffuses in the opposite direction, counter to the direction of its own complementary concentration gradient.Note also that the cubic source term has also appeared in the context of population genetics, as a correct modification to Fisher's equation for gene frequencies when the genes are not Mendelian but are neither fully dominant nor fully recessive [12][13][14].
The usual simplified Cahn-Hilliard equation is: The operator ∇ 2 is the Laplacian and ∇ 4 is the biharmonic operator.The negative Laplacian term (with b < 0) is destabilising, and generates time-reversed diffusion.It has positive spectrum, with amplification of perturbations at high spatial frequency.However, when the spatial frequency is high enough, the instability is restrained by an opposing negative biharmonic term (with a < 0) which has a negative spectrum, and at high spatial frequencies, the overall spectral values are negative rather than positive.
Although these diffusion equations do not minimise any action functional, they do result ultimately from an energy functional.In the Cahn-Hilliard approach, the phase flux is the negative gradient of a chemical potential multiplied by a conductivity (e.g., [4,6,9,15,16]).The chemical potential Φ is the variational derivative Φ = δE /δθ, where E is the energy functional.Further, it is reasonable to include within E a component of strain energy B(θ)|∇θ| 2 that depends locally on the phase field value as well as on its squared gradient.Karali and Katsoulakis [16] showed that in some cases, a combined Allen-Cahn-Hilliard equation could indeed be derived from an energy principle, allowing for both a gradient flow and a reaction term that allows for relaxation towards equilibrium: Comparatively few exact solutions are known for time-dependent multi-dimensional fourth-order nonlinear reaction-diffusion equations.Galaktionov and Svirshchevskii (Example 6.72 of [17]) constructed some solutions for the thin film equation with the absorption term −W (θ) = − (θ), allowing for θ simply to be quadratic in Cartesian coordinates x i .These solutions, along with a number of others, were derived by Cherniha and Myroniuk [18] in a comprehensive Lie symmetry classification.In the current article, some more realistic explicit solutions with meaningful boundary conditions will be constructed after a reduction by a strictly nonclassical symmetry.The solutions have nontrivial spatial structure but they approach the pure phase θ = 0, exponentially in time.
In Section 2, we show that a nonclassical symmetry allows one to separate variables to linear equations in space and time, when the nonlinear reaction and diffusivity obey a single relationship.This allows us, in Section 3, to construct solutions that approach the phase θ = 0 asymptotically in time when the Allen-Cahn source term is exactly the cubic function that was given in [5].
In Section 4, we show how these fourth-order reaction-diffusion equations relate to an energy principle for a flux potential variable that is the analogue of the matrix flux potential that is well-known from nonlinear porous media studies [19,20].
In Section 5, we conclude with a discussion of the solutions, and some open problems.

Role of a Fourth-Order Kirchhoff-Helmholtz Equation
For the purpose of reducing the number of variables and producing special symmetric solutions, the requirement of leaving the whole solution manifold of a PDE invariant could be weakened so that only a subset of solutions is invariant under partial or conditional symmetry [21].In particular, if a symmetric solution can be constructed, then by definition, the target equation must be compatible with the invariant surface condition that is simply an expression of that invariance.This led Bluman and Cole [22] to the definition of nonclassical symmetry as a one-parameter transformation, an analytic function of a single real parameter connected to the identity transformation at = 0 that leaves invariant the system consisting of a target PDE together with the invariant surface condition.The first example to be considered [23] was the linear heat equation supplemented by the invariant surface condition: Subscripts will denote partial derivatives with respect to the subscripted variables.A full set of determining equations was constructed for the coefficients (ξ 1 , ξ 2 , η) of the infinitesimal symmetry generator: Unlike the determining relations of the classical symmetries that are linear PDEs, the nonclassical determining relations for (ξ 1 , ξ 2 , η) are nonlinear.The nonclassical determining relations for the linear heat equation were not solved for three decades [24].One class of nonlinear PDEs that has been fully classified by solving the nonclassical determining relations is the class of semi-linear (1 + 1)-dimensional reaction-diffusion equations with a linear diffusion term and general nonlinear reaction term [25,26].A more extensive account of non-classical symmetry classification is given in the recent book by Cherniha and Davydovych [27].
Pertinent to the current study, a general nonclassical symmetry classification of fully nonlinear reaction-diffusion equations in (2 + 1)-dimensions was completed in [28].The class of target PDEs was: It is convenient to change the dependent variable to the Kirchhoff [29] variable: so that the class of PDE ( 4) is expressed as Not only does this reduce the number of terms in the PDE but it also simplifies the flux to J = −∇u.Since θ = 0 is a stable stationary state, θ 0 is henceforth chosen to be 0.
Without loss of generality, nonclassical symmetries have ξ 0 = 1 or ξ 0 = 0. Hence, in [28], the system to be analysed was (6), combined with: Among the compatible nonclassical symmetry generators, there was the unexpected simple combination: Note that if D(θ) is specified, then (9) gives R(θ) explicitly.However, if R(θ) is specified, then (9) must be treated as a differential equation to solve for u(θ), after which D(θ) = u (θ).The symmetry reduction gives: Provided the nonlinear diffusivity function and the nonlinear reaction function are related as in ( 9), the nonlinear reaction-diffusion equation is amenable to separation of variables to two linear equations, allowing an arbitrary solution of the Helmholtz equation (Equation (10)).This led to the only known exact time-dependent solution of Arrhenius combustion with diffusion [30] by a symmetry reduction that applies in n ≥ 1 spatial dimensions.
Until recently, it was not apparent that the non-classical reduction could be applied to fourth-order reaction-diffusion equations.At this point, it is important to note that provided (9) holds, the same nonclassical symmetry applies when the Laplacian operator in the PDE ( 6) is replaced by any linear operator, that is, Lu can be any time-independent linear combination of u and its spatial derivatives, In a recent article [31], this was applied directly to a soil-water-crop system, by choosing L to be the Kirchhoff operator Lu = ∇ 2 u − u x 3 .Now we choose Lu = a∇ 4 u + b∇ 2 u with a, b ∈ R, so that u and θ satisfy: A direct substitution of: into (12) leads to the single requirement (9).Equation ( 14) could be regarded as a fourth-order Helmholtz equation.
The constrained system consisting of (6) with R(θ) and D(θ) satisfying ( 9), supplemented by the simple invariant surface condition u t = Au, has an infinite dimensional point symmetry group: where Φ(x, t) is an arbitrary solution of the linear Equation ( 14).However, this type of symmetry does not appear in the classical symmetry group of the original unconstrained PDE (6) that is not integrable.
In the sense that we can construct an infinite dimensional linear space of solutions of a restricted form, the unconstrained PDE is semi-integrable.The linear space of constructed solutions originates from a linear PDE in two independent variables compared to the original three independent variables.

Amenable Diffusivity and Reaction Functions
The amenable diffusivity function for the second-order multi-dimensional nonlinear reaction-diffusion equation with the cubic reaction term was already investigated in [32] in the context of population genetics.However, in that application the steady state θ = 0 (meaning no advantageous genes present) is unstable, whereas in the current application, θ = 0 is a stable pure phase.
First consider the canonical cubic reaction term that appears in (1), We assume that the phase that is approached exponentially is θ = 0, requiring A = −|A| < 0. In a phase pattern formation, there is some interest in the case that the fourth-order diffusion is dissipative (aD < 0) whereas the second-order diffusion is backwards in time (bD < 0), giving a negative contribution to entropy production [33].We will focus on this case below, and assume that these sign properties are determined by the parameters a and b taking negative values after assigning D > 0.
At this point it is convenient to introduce dimensionless variables.We assume D(θ) has the usual dimensions of diffusivity, This inverse length is around the critical wave number, above which the fourth-order diffusion overcomes the second-order backward diffusion, resulting in attenuation and stability.A natural choice of length and time scales is s = √ a/b and t s = 1/s.Define dimensionless coordinates (x * , y * , z * ) = (x, y, z)/ s and t * = t/t s , with dimensionless Laplacian ∆ * = 2 s ∆.In that system of units, s is scaled to 1, whereas the time scale for the particular solution to decrease exponentially is written as 1/|A * | with A * = At s .Defining: Equation ( 13) reduces to: With A * = A/s and κ * = κ|a|/b 2 , Equation (9) can also be written in using dimensionless variables, without a, b, or s occurring explicitly, as can other equations of interest above.From here on, non-dimensional coordinates will be assumed, and asterisks will be dropped.
Since R(0) = 0 and u(θ) = 0 + D(0)θ + O(θ 2 ), (9) implies: In summary, From here on, we require |A| > 1/2.Naturally, the rate of exponential decay (with time scale 1/|A|), will be related to the strength of the sink term (with time scale 1/s).Equation (19) shows that the diffusivity must be increasing from θ = 0 to some value θ m ∈ (0.5, 1) where D reaches a maximum.It is physically reasonable to expect that the mobility of one phase through a mixture will change, resulting in a change of diffusivity, as the composition of the mixture changes.For example, take |A| = 5/2.This results in D(0) being 4|A|/(5κ), which is 20% lower than D(0.5).
The required solution of differential Equation ( 9) is the fixed point of the map: 2 )/κ (constant).Then: This already has the correct values at the three locations given in (19), however we now require |A| > 9/16 to avoid singularities in the region θ ∈ (0.5, 1).From D 1 (θ) we can calculate the corresponding: (23) so that the next iterate D 2 (θ) can also be calculated explicitly from (20).
As demonstrated in [32], these iterates may converge rapidly to the exact solution.In this example, although the analytic expressions for D 1 (θ) and D 2 (θ) are quite different, numerically they may be practically indistinguishable, as illustrated in Figure 1 Given the cubic reaction term ( 16), D 1 (θ) is already a good approximation for the diffusivity function that leads to our exact solution for the phase field.Conversely, we may choose the diffusivity function to be D 1 (θ) and then exactly construct a reaction function R 1 (θ) that has the same generic properties as the canonical cubic function, but which exactly leads to the analytic solution for the phase field.From (9), Figure 2 plots the cubic R(θ) along with R 1 (θ) for A = −5/2.

Interior Solutions for Slabs, Cylinders, and Spheres
In order to gain some understanding of the nature of the solutions, we consider interior solutions for slabs, cylinders, and spheres in dimensions n = 1, 2, and 3 respectively.The four independent solutions will be trigonometric and hyperbolic functions (n = 1), Bessel and modified Bessel (n = 2), and spherical Bessel and modified spherical Bessel functions (n = 3).The solutions decrease uniformly in t towards θ = 0.Even if the initial condition has θ > 0.5 in some significant sub-domain, the region where the source term is tends to stabilise phase θ = 1.
Therefore, this global approach to Phase 0 should be controlled by the boundary conditions which may include, for example, ideal contact with the exterior pure phase at θ = 0, or a Robin boundary condition for extraction of Phase 1.
In general coordinate systems with our adopted dimensionless variables, the simplified form of Equation ( 14) factorizes as: where: For radial flow within the ball ∑ n i=1 x 2 i = r ∈ (0, r 0 ), the Laplacian operator within ( 25) is . The two operator factors commute so (25) reduces simply to a system of two alternative second-order Helmholtz equations.The general solution Φ(r) of ( 25) in one, two, and three dimensions, with four free parameters c i ∈ R, is: Here, the spherical Bessel functions are defined simply as j 0 (r) = sin(r)/r, y 0 (r) = cos(r)/r, i 0 (r) = sinh(r)/r and k 0 (r) = cosh(r)/r.
For radial flows, the flux density of Phase 1 is: For all values of n we impose zero flux at r = 0.For n = 3, considering the asymptotic behaviour of J(r) as r → 0 shows that J(0) = 0 imposes two conditions on c 2 and c 4 : These admit no non-trivial solutions, so c 2 = c 4 = 0.For n = 3 alternatively imposing either Φ (0) = 0 or Φ (0) = 0 also gives two conditions on c 2 and c 4 similar to those above, which also imply c 2 = c 4 = 0.
For n = 1, J(0) = 0 only implies the condition: As Φ(0) is finite, the n = 1 case does admit physical solutions where J(0) = 0 and Φ (0) = 0. We can eliminate this extra degree of freedom by imposing Φ (0) = 0, to match the admissible solutions discussed above for n = 2 and n = 3.This gives an extra condition on c 2 and c 4 : so that again only c 2 = c 4 = 0 is possible.The boundary r = r 0 is in contact with pure Phase 0. Imposing Φ(r 0 ) = 0, implies: We have found in examples that the zero-flux boundary condition at r = r 0 will result in an inadmissible solution that is not positive semi-definite.However there is a solution that connects smoothly with an external slab of Phase 0, beginning at the boundary r = r 0 .Assuming Φ (r 0 ) = 0, we deduce: This has a sequence of solutions for r 0 .With κ = 1, the first of these are r 0 ≈ 2.0616, 2.7012 and 3.2653 for n = 1, 2, and 3 respectively, showing an significant increase with dimensionality.Higher values in the sequence of solutions for r 0 are found to produce unphysical solutions exhibiting negative values of u.With diffusivity D 1 (θ) and reaction term R 1 (θ) shown in Figures 1 and 2 respectively, θ is given explicitly in terms of u as: Although solutions for u(x, t) at different times are geometrically similar, this is not true for θ(x, t) which is related to u by a nonlinear transformation.The solution for θ(x, t) is depicted in Figure 3 for A = −5/2, κ = 1.The solutions in dimensions n = 1, 2, 3 are close after the three different domain radii r 0 are rescaled to the same value.The solutions depict an initial region in the neighbourhood of the origin, that consists predominantly of Phase 1, but is connected continuously through a thin mushy zone to a surrounding block of Phase 0, at r = r 0 .This mushy region is small enough to be stable under the Cahn-Hilliard evolution.Under the influence of the reaction, Phase 1 is normally stable, but its decay in this case is driven by the boundary conditions at r = r 0 , a boundary where the phase mixture is connected smoothly to an exterior bank of Phase 0. For example, from t = 0 to t = ∞, the exact radial solution on the interior of the sphere, depicted in Figure 3, has 90% of the initial volume of Phase 1 material exiting through the boundary, while the remainder undergoes a change to Phase 0 within the interior.
Unphysical solutions, with negative values of u, were found above for large values of r 0 .If one attempts to find the largest possible value of r 0 that produces a non-negative solution for a given κ, it soon becomes clear that for the class of solutions considered above, this corresponds to the first non-zero solution r 0 , of the above conditions Φ (r 0 ) = 0 shown in (33).This is a consequence of our solutions being a linear combination of one oscillatory, and one monotonically increasing function that take finite values at r = 0.As such, Figure 3 shows scaled examples of solutions with the maximum possible value of r 0 .Selecting smaller values of r 0 produces solutions with a non-zero negative slope at r = r 0 .Figure 4 shows variation of the maximum possible r 0 value as κ varies.In all cases max r 0 → 0 as κ → ∞.The case κ = 0 does not produce valid solutions via our methods above, however the largest possible values of r 0 are obtained as κ → 0, and Equation (33) approaches the following simple conditions in one, two, and three dimensions respectively:

Interior Solutions for Rectangular Domains
Under idealised controlled boundary conditions on simple slab, cylindrical, or spherical domains, 2r 0 gives a maximum achievable diameter of the transition zone.However, there are many other solutions with finer sub-structure.This is consistent with the Cahn-Hilliard critical wavelength being of order 1, in dimensionless terms.The initial condition for u(r, 0) = φ(r) may be any positive solution of the linear fourth order Helmholtz equation.
If we adopt Cartesian coordinates and restrict our attention to solutions that are identically zero on the boundaries |x| = x 0 , |y| = y 0 , the dimensionless form of Equation ( 14) with a, b < 0 can be solved via the separation of variables.Setting φ(r) = X(x)Y(y) we find two families of solutions where either: If we focus on the latter case, and impose Y(±y 0 ) = 0 we find compatible values of λ: Meanwhile X(x) must satisfy the differential equation: provided that λ m < 1 2 + κ + 1 4 .This produces a solution similar to the first of ( 27): The second form of X m (x) above has the boundary condition X m (±x 0 ) = 0 imposed.Note that there is a maximum possible value of x 0 that will produce a non-negative X 1 (x) function-hence, the situation is qualitatively similar to the examples of φ(r) in various dimensions discussed in Section 3.
On the other hand, when λ m > 1 2 + κ + 1 4 we can define: The corresponding solution X m (x) then loses all periodicity: The above form of X m (x) will apply for all m = 1, 2, . . .if: that is, if y 0 is sufficiently small.When this condition is met there is no restriction on the size of x 0 implied by needing a non-negative solution X 1 (x).
In summary, if y 0 is sufficiently small we can choose x 0 to be as large as we want, however if y 0 > y 0,c shown above, there is some upper limit to the size of x 0 that will produce valid solutions.When all separable modes are accounted for, we can calculate the boundary separating feasible combinations of x 0 and y 0 from infeasible combinations.Figure 5 illustrates this boundary for various values of κ.Note that as y 0 → ∞ the critical values of x 0 approach those indicated by the n = 1 curve of Figure 4.When considering combinations of x 0 and y 0 that admit physical solutions, the eigenvalue corresponding to m = 1 with c 6 = 0 is of primary importance as both X 1 (x) and Y 1 (y) are non-negative.We view solutions with finer substructure as alterations of this solution that maintain the non-negativity of φ(x, y).
An example with a pronounced sub-structure is shown in Figure 6.Here, κ = 1, x 0 = 0.3, y 0 = 1.2, and φ(x, y) is proportional to X 1 (x)Y 1 (y) + 0.225X 5 (x)Y 5 (y) where: From a brief exploration of the parameter space, it appears that solution regions with aspect ratios close to 1 typically allow little variation from X 1 (x)Y 1 (y), and tend to be unimodal.Solution regions that are more elongated can be multi-modal, and allow noticeable sub-structure as in the example of Figure 6.We can also consider rectangular solutions that are unbounded in a single direction, similar to the one-dimensional solution of Section 3. Let us begin by assuming that boundary conditions X(±x 0 ) = 0 are applied to Equation (38) for X(x).The eigenvalues of Y (y) + λY(y) = 0 are still restricted by requiring a non-negative solution that remains finite as y → ±∞.Hence the separable solution must include k 0 = λ 0 = 0, Y 0 (y) = 1 as its foundation.Equation (40) then gives a non-negative X 0 (x), provided that x 0 is not larger than the critical values illustrated in Figure 4 for n = 1.For values of x 0 significantly less than these critical values, we can produce noticeable sub-structure by adding other modes to this base solution while ensuring that φ(x, y) remains non-negative and finite.Figure 7 shows one such solution with κ = 1 and x 0 = 0.4.It is proportional to X 0 (x) − 0.3X λ (x) cos( √ λy), with λ = 25/16.Here X 0 (x) is given by (40) with c 5 = 1, c 6 = 0, and ω m and λ m replaced by ω and λ as given by Equation (26), respectively.X λ (x) is also given by (40) with c 5 = 1, c 6 = 0, but with ω m and λ m replaced by: It is also possible to obtain a rectangular unbounded solution by first applying Y(±y 0 ) = 0 on the solutions of Y (y) + λY(y) = 0. Seeking a finite non-negative value of X(x) again supports the conclusion that no solutions are possible for y 0 values greater than those indicated by the n = 1 curve of Figure 4.
In summary, our consideration of solutions φ(x, y) interior to a rectangular region is compatible with the maximal values of r 0 derived in Section 3, and shows that exact solutions with varied substructures are possible.

Energy Formulation
Let us consider the original equation in (13), Equation ( 46).In more physically relevant situations R represents a competing energy contribution usually in the form of a non-negative multi-well function in the variable u.When −R is a non-convex non-negative polynomial it may still be possible to ensure the existence (although not the uniqueness) of the solution u to the minimum problem for I(u) possibly in a subset of H 2 ϕ (Ω) and characterize u as a solution to Equation (48).Instead, when b < 0 and −R = −cu 2 with c > 0 we have genuine competition in the summands of the energy I and in turn I may fail to be positive-definite.Under these circumstances, solutions to Equation (48) correspond to various critical points of I including saddles or more intricate situations.In those cases, existence theorems should be investigated by methods possibly based on min-max principles.
An alternative approach (e.g., [15,16]) is to express the phase flux density as being proportional to the negative gradient of a chemical potential energy density that is itself the variational derivative of a total internal energy functional.This comes about because the drift velocity of the mobile particles of Phase 1 quickly reaches a terminal velocity that enables the driving gradient of chemical potential energy to be balanced by a mechanical resistance force that is proportional to velocity and in the opposite direction.The new element that we need here to justify the non-standard Equation ( 6) is that the variational derivative is taken with respect to the flux potential u, rather than the phase variable θ.Then we choose: In this set-up, R(θ) must be regarded as an externally driven Phase 1 production term that assists relaxation towards the stable fixed points.

Conclusions
A nonclassical symmetry reduction has led to some exact positive solutions for multi-dimensional nonlinear reaction-diffusion equations with both fourth-order diffusion and second-order backward diffusion.In all likelihood, these are the only known exact solutions for equations of the Allen-Cahn-Hilliard type that have meaningful boundary conditions.The presented solutions have a zero flux and a zero gradient at one central boundary point, plus a prescribed value, representing a pure phase, as well asa zero gradient at the other boundary.The initial condition has the phase field taking the value θ = 1 at r = 0.Even though Phase 1 is stabilised by the reaction term, the boundary conditions force the solution to have the system uniformly approaching Phase 0. An interior cylindrical or spherical inclusion with a mixed phase could not withstand the inward advance of Phase 0 from the outer boundary.
Structured phase domain patterns have been seen in numerical solutions of the Cahn-Hilliard equation (e.g., [35]).Some exact solutions with pronounced sub-structure have been produced here too.Since the exact solution involves a general solution of the fourth-order linear Helmholtz equation for the flux potential, it is possible that more complex patterns could be incorporated in the future.Any solution of the fourth-order Helmholtz equation can be used in this method, even those that apply in more complicated asymmetric domains.
We remark the approach based on the Helmholtz substitution adds additional structure to the model in terms of a variational principle as the by-product of symmetry and linearity.Although the thorough investigation of the relation between the non-linear model in the unknown variable θ and the transformed model in u is not the scope of the present contribution, in Section 5 some insight is provided for some special situations.The construction used here became possible because of a nonclassical symmetry that applies when the nonlinear diffusivity and the nonlinear reaction term together obey a single relationship.Then an additional constraint θ t /D(θ) = A θ θ 0 D(s)ds, added to the governing PDE, results in an integrable system.Reduction of variables then results in the linear fourth-order Helmholtz equation.In that sense, the original multidimensional nonlinear scalar PDE is partially integrable because one additional constraint leads to a linear PDE with one fewer independent variables.It remains to be discovered whether nonclassical symmetry classification can unearth other semi-integrable PDEs-possibly defined by different fourth-order non-linear operators-in this way, and more generally, how compatible constraints may be selected.

Figure 5 .
Figure 5. Critical combinations of x 0 and y 0 separating the feasible rectangular solution region from the region where x 0 and y 0 in combination do not admit physical solutions.