## Abstract

We consider a simplified model of a two-phase elastic solid with small interfacial energy in forced anti-plane shear. We propose a novel approach to finding meta-stable equilibria (local minima of the potential energy), based upon global bifurcation theory, *a priori* bounds and numerical path following. In particular, we treat the reciprocal of the capillarity coefficient as a continuation parameter, and we fully exploit the hidden symmetry of our elliptic boundary-value problem, both strategies of which are crucial for obtaining rigorous existence results and for performing reliable and efficient computation. We obtain branches of locally stable equilibria exhibiting phase nucleation (and anti-nucleation) and fine layering of mixtures.

## 1. Introduction

In recent years, rigorous mathematical tools have emerged for the prediction of microstructure for a broad class of shape-memory alloys. Within the context of sharp-interface models, the philosophy there is one of global energy minimization—or rather the construction of weakly convergent minimizing sequences (e.g. Ball & James 1987, 1992; Bhattacharya 2003). For materials characterized by a free energy having multiple wells of equal height and in the absence of loading, the approach predicts such things as the inclinations of twinning planes and martensite–austenite interfaces and the correct volume-fraction of phase mixtures with striking success. However, sharp-interface models suffer from the defect of allowing both infinite refinement and rearrangement of microstructure. Also, the above-mentioned approach typically ignores equilibrium conditions (across sharp interfaces), does not easily generalize to problems in the presence of loadings and/or energy wells of unequal height, and does not directly address local minima of the energy. The study of energy minimizers of two-phase problems in the presence of asymptotically small interfacial energy or capillarity is difficult, but there are results in mostly one-dimensional and special two-dimensional settings pertinent to solids (e.g. Carr *et al*. 1984; Müller 1993; Kohn & Müller 1994; Huo & Müller 2003; Yip 2006).

Experiments on multi-phase solids are often carried out by slowly loading an initially homogeneous specimen while carefully observing the emergence of increasingly exotic phase mixtures, patterns or microstructure (e.g. Chu 1993). Methods of global bifurcation/continuation theory and symmetry-breaking techniques are well suited for the analysis of such problems. The idea is to obtain connected branches of *equilibria* and determine their stability as one or more loading parameters evolve. With the exception of various one-dimensional models (e.g. Lifshitz & Rybnikov 1985; Triantafyllidis & Bardenhagen 1993; Vainchtein *et al*. 1999), and also some studies of the two-dimensional Cahn–Hilliard problem (Kielhöfer 1997; Maier-Pappe & Miller 2002), bifurcation/continuation methods have not had much of an impact on the field. Two possible reasons for this are: (i) realistic models incorporating small interfacial effects reveal localization phenomena, which are typically not captured by well-known asymptotic bifurcation techniques. More specifically, in the vicinity of a bifurcation point, the non-trivial solution behaves much like the eigenfunction of the linearization. Indeed, the above-mentioned studies using global bifurcation methods all demonstrate that the relevant physical solutions are global, i.e. ‘far’ from the trivial solution; (ii) the global bifurcation analysis of the higher-order pde's and systems of pde's appropriate for two- and three-dimensional models is much harder than that of ode's and second-order pde's.

We consider a two-dimensional problem of the former type in this paper, *viz*. a simplified anti-plane shear model of two-well elasticity in the presence of isotropic interfacial energy with hard loading on the boundary. In the absence of interfacial effects and loading, this model was first introduced in Rosakis (1992) for the study of meta-stable equilibria of crystalline solids with cubic symmetry. Other models for two-phase solids in anti-plane shear have been proposed and studied in Fried (1993) and Pettinger & Abeyaratne (2000). Here, we are led to the study of a fourth-order elliptic pde over a two-dimensional domain. In particular, the model is more complex than that of the steady-state Cahn–Hilliard model, due to the fact that the two-well potential is a function of the strain *vector* field.

As in Maier-Paape & Miller (2002), we first divide the governing equilibrium equation by the small capillarity coefficient (appearing in front of the principle part of the operator), treating its reciprocal as a bifurcation/continuation parameter, in addition to the natural hard-loading parameter representing the magnitude of the externally applied deformation on the boundary. We first demonstrate and exploit hidden symmetry in our boundary-value problem. This enables a global symmetry-breaking bifurcation analysis from the homogenous trivial solution, with each branch of solutions having precisely enumerated symmetry properties. Next we establish *a priori* bounds on all possible solutions. When combined with the global bifurcation results, this yields the rigorous existence of solutions characterized by small interfacial effects. In particular, unbounded branches necessarily contain solutions for arbitrarily small capillarity, the latter of which are excellent candidates for meta-stability. We then perform numerical path following, with the reciprocal of the capillarity coefficient as the path parameter, while seeking local minima of the potential energy. Once found, we then fix the value of capillarity coefficient and continue the path following in the actual loading parameter. In this way, we compute branches of locally stable equilibria exhibiting phase nucleation (and anti-nucleation) and fine layering of phase mixtures—all in qualitative agreement with experiment (e.g. Chu 1993).

The outline of the paper is as follows: in §2, we first formulate the problem on a square domain. We then reformulate the pde as an O(2) equivariant bifurcation problem and then perform a global bifurcation analysis yielding the existence of solutions branches in the large. In §3, we obtain *a priori* bounds on all solutions, resulting from a combination of positivity, energy estimates, embedding theorems and elliptic regularity. In §4, we discuss our numerical path-following strategy. Here, we fully exploit the symmetries of each solution by computing on reduced domains. Of course, the check for local stability is carried out for the extended solution on the full domain.

## 2. Global bifurcation analysis

We consider an infinite cylindrical body of square cross-section, , undergoing anti-plane-shear displacement, , in response to imposed boundary displacements of the form *w*|_{∂Ω}=*λx*, where *λ* is a loading parameter. The total potential energy of our model is the following:(2.1)where ∇*w* and ∇^{2}*w*≔∇(∇*w*) denote the first and second displacement gradients, respectively. Also, is a non-convex function with two equal-energy wells, corresponding to different phases or twin variants. The small parameter *ϵ*>0 is the capillarity coefficient of the regularizing, second-gradient term, the later of which in equation (2.1) models interfacial energy. More specifically, we choose(2.2)The associated Euler–Lagrange equilibrium equation associated with equation (2.1) is(2.3)where *u*≔*w*−*λx*; *κ*≔1/*ϵ*; *f*(*t*)≔6*t*(1−*t*), Δ(.) denotes the Laplace operator; Δ*u*=*u*_{xx}+*u*_{yy} and denotes the biharmonic operator. Observe that equation (2.3) admits the trivial solution *u*≡0 for all .

We wish to establish the existence of global solution branches bifurcating from the trivial solution. Obviously, we are interested in solutions for which *κ* is large. As a first step in the bifurcation analysis, we formally linearize (2.3) about *u*≡0,(2.4)which admits the non-trivial solutions(2.5)provided that(2.6)

While equations (2.5) and (2.6) suggest the possibility of bifurcation from the homogeneous solution *w*≡*λx* (*u*≡0), the lack of smoothness of the domain *Ω* precludes a straightforward rigorous analysis. We can easily get around this by observing and exploiting the hidden symmetry inherent in our problem. We first consider the partial differential equation in (2.3) on the infinite strip , which is readily shown to admit the following equivariant symmetries:(2.7)Following the approach in Healey & Kielhöfer (1991), it is routine to demonstrate that for each , the pde in equation (2.3) defines a mapping(2.8)where the function spaces are defined as follows: for *k*≥2,(2.9)where *C*^{k,α}(*S*), *k*=0, 1, 2, …, denotes the usual space of *k*-times (locally) Hölder continuously differentiable functions. Clearly, for each , and are Banach spaces when equipped with the appropriate uniform Hölder norms over the closure of the subdomain . Finally, observe that any solution of(2.10)with *F* restricted to as in equation (2.8), yields a solution of equation (2.3). Similarly, any solution of the boundary-value problem (2.3) belonging to can be used to generate a solution in (via equation (2.7)_{2} and 2*π*-periodic extension). Accordingly, we henceforth analyse equation (2.10) on the periodic Hölder spaces as in equation (2.8). In particular, the linearization (2.4), now corresponding to the total derivative of equation (2.10), along with equations (2.5) and (2.6) are now rigorous.

Although it is straightforward to perform a two-parameter global bifurcation analysis of equation (2.10), we instead emphasize two one-parameter choices in anticipation of our forthcoming numerical results. However, we first point out our interest in thin ‘stripes’, typically associated with twinning (e.g. Ball & James 1987; Chu 1993). In particular, we henceforth focus exclusively on the case *m*=1 in equations (2.5) and (2.6). For any fixed value *κ*=*κ*_{0}, the characteristic equation (2.6) yields(2.11)Observe that for a given , equation (2.11) yields a pair of critical values that are real only for *κ*_{0} sufficiently large (*ϵ*_{0}=1/*κ*_{0} sufficiently small), which is precisely the same situation that occurs in the analysis of one-dimensional and two-dimensional Cahn–Hilliard models (cf. Lifshitz & Rybnikov 1985; Triantafyllidis & Bardenhagen 1993; Kielhöfer 1997; Vainchtein *et al*. 1999; Maier-Paape & Miller 2002). Specifically, the smaller the value of *ϵ*_{0}, the more ‘modes’ (2.5) participate. Alternatively, we may fix *λ* and solve for critical values of *κ*. Here we choose *λ*=1/2 (at the peak of the so-called spinodal region), in which case equation (2.6) leads to(2.12)

It is easy to show that each of the points and (*κ*_{n}, 0) correspond to simple bifurcation points in the fixed-point spaces . Moreover, each of the local bifurcating solution branches are ‘pitchforks’, and each solution is asymptotic to the appropriate eigenfunction (2.5) (with *m*=1). Each of the pitchforks corresponding to the points and (*κ*_{n}, 0) are ‘supercritical’, while those corresponding to the points are subcritical. Throughout this work, we adopt the minimum energy criterion of stability, i.e. an equilibrium solution of equation (2.3) is *stable* if it renders the potential energy (2.1) a local minimum. Otherwise, the solution is said to be unstable. Although we do not give the details, it can be shown that *all* of these non-trivial local solutions are unstable (failing the second-variation condition) with the exception of the branch corresponding to (*κ*_{2}, 0). However, even that (local) solution is also physically irrelevant—due to the relatively small value of in a neighbourhood of the bifurcation point.

Accordingly, we next perform a (standard) global bifurcation analysis, with the eventual hope of finding stable equilibria. Observe that the linear operator , defined by for all , has a bounded inverse denoted by *L*^{−1}. Moreover, is compact. Hence, equation (2.10) is equivalent to the operator equation(2.13)where is the compact map defined by(2.14)A direct application of a well known theorem of Rabinowitz (1971) now yields:

*Each pair* , *with κ*_{n} *given by equation* (2.12), *is a bifurcation point of a global continuum of non-trivial solution pairs of equation* (2.13) (*with λ*=1/2), *denoted by* , *n*=2, 3, …. *Each branch Σ*_{n}, *n*=2, 3, …, *is either unbounded in* *and*/*or meets the trivial line of solutions at some other point* , *with κ*_{*}≠*κ*_{n}. *Likewise*, *each of the pairs* , *assuming that* *given by equation* (2.11) *are real*, *is a bifurcation point of a global continuum of non-trivial solution pairs of equation* (2.13) (*with κ*=*κ*_{0} *fixed*), *denoted by* , *each of which is characterized by the two* ‘*Rabinowitz alternatives*’, *as described for Σ*_{n}.

## 3. *A priori* bounds

In this section, we establish bounds on all possible solutions of the partial differential equation (2.3) as a mapping on the periodic Hölder space , in terms of the Hölder norms, denoted by , for *k*=2 and 4 (observe that *Ω*_{1}=*Ω*), cf. theorem 3.4. Along the way, we also make use of the usual *L*^{2}(*Ω*) inner-product, , and the *L*^{p}(*Ω*) and *W*^{j,p}(*Ω*) norms, denoted and , respectively, *p*>1, *j*=1, 2, …. Generic constants are denoted throughout by, e.g. *C*_{λ}, the subscript of which reflects dependence upon *λ*. All constants typically depend upon the domain *Ω* as well, but not upon the solution *u*. Our first preliminary result is:

*There exists a positive constant C*_{κ,λ} *such that*(3.1)

We multiply the pde in equation (2.3) by *u* and integrate over *Ω* yielding(3.2)where *F*(*t*)=3*t*^{2}−2*t*^{3} is the anti-derivative of *f*(*t*). (Of course this is simply (2.1) again—written in terms of equation (2.2) and *u*≔*w*−*λy*.) Observe that the left side of equation (3.2) is bounded below by 〈*u*_{x}, *u*_{x}〉, and note that the boundary conditions in equation (2.3) imply that . These observations combined with equation (3.2) yield(3.3)By the Cauchy–Schwarz inequality, equation (3.3) leads to(3.4)where we have used the fact that the minimum value of 6*λ*^{2}−6*λ*+1, , is −1/2. From this we conclude(3.5)But then (since *F* is a cubic polynomial) we see thatand by virtue of equation (3.2) we have(3.6)From equation (3.6), it follows that the *L*^{2} norms of the gradient and second gradient of *u* are each bounded above by a constant—for the latter, the constant depends upon both *κ* and *λ*. Using Poincare's inequality, we then deduce(3.7)In particular,(3.8)and by embedding (Adams 1975) we deduce(3.9)for all *q*≥2. We choose *q*=12 to obtain (*f* is a quadratic polynomial)(3.10)Finally, Hölder's inequality (3.7) and (3.10) now yield

Some easy consequences of arguments in the proof lemma 3.1 are:

*For κ*=0, *u*≡0 *is the only classical solution of problem* *(2.3)*. *Similarly, u*≡0 *is the only classical solution of equation* *(2.3)* *when λ*=−1/6 *and when λ*=7/6 *for all κ*>0.

When *κ*=0, equation (3.2) implies that ∇^{2}*u*≡0, and the first claim then follows from the boundary conditions. Next, assume that *u* is a non-trivial solution of equation (2.3) for *λ*=−1/6. From equation (3.3), we then have(3.11)Clearly, the left side of equation (3.11) is strictly bounded from below as follows:

On the other hand, the arithmetic–geometric means inequality applied to the right side of equation (3.11) yieldsThus, equation (3.11) leads to a contradiction unless *u*≡0. The same argument goes through when *λ*=7/6.

For our next result, we need the following elliptic *L*^{p} regularity result due to Agmon (1960, theorems 7.1' and 8.1):

*Let* *be a bounded domain with a smooth boundary* ∂*Θ*, *and let* *denote a uniformly strongly elliptic operator*, *with* . *Let A*^{*} *denote the adjoint operator of A*, *assume that g*∈*W*^{i,p}(*Θ*) *for some fixed integer i*∈{0, 1, …, *j*}, *and suppose that v*∈*L*^{q}(*Θ*) (*q*>1) *satisfies* (*the weak form*)(3.12)*for all* (*the set of all infinitely differentiable functions on Θ of compact support*), *where* 〈., .〉_{Θ} *denotes the L*^{2} *inner-product over the domain Θ*. *Then*,(3.13)*where the constant C is independent of u and g* (*but dependent upon the coefficient functions comprising A and the domain Θ*).

Actually, proposition 3.1 is not directly applicable to our problem, ostensibly due to the fact that the square domain *Ω* has corners. However, the proof of theorem 8.1 in Agmon (1960) relies upon the ability to homeomorphically map small half-open neighbourhoods of points on the smooth boundary into hemispheres. Excluding small neighbourhoods of the four corners, we can do this directly with ‘half-discs’ along ∂*Ω* as required in Agmon (1960). For the corners, we use the symmetry inherent in equations (2.9) and (2.10) as follows. We cover each corner with a sufficiently small ‘quarter-disc’ neighbourhood, denoted by *Q*_{j}, *j*=1, …, 4. We let *H*_{j} denoted the half-disc obtained by reflecting *Q*_{j} about the appropriate axis (either the *y*-axis or the line *x*=*π*). Then, the periodicity and oddness properties of immediately yieldensuring that the same argument leading to inequality (8.4) of Agmon (1960) holds. In particular, proposition 3.1 applies to our periodic problem in terms of the parent domain *Ω*.

First, observe that equation (2.3) can be rewritten as(3.14)for all , where and . Then equation (3.13) of proposition 3.1 with *i*=0 and *p*=3/2 yieldsand from equations (3.1) and (3.7) we then conclude(3.15)Our main result of this section is the following:

*For any solution u of equation* *(2.10)* *there are positive constants C*_{λ,κ} *and K*_{λ,κ}, *depending on λ*, *κ and Ω*, *but independent of u*, *such that*(3.16)*and*(3.17)*for* 0<*α*≤2/3.

Observe that *u* trivially satisfies the equation(3.18)for all . Inequality (3.13) of proposition 3.1 with *i*=2 and *p*=3/2 now delivers(3.19)where the second inequality follows from equations (3.7) and (3.15). The claim (3.16) now follows from a well-known imbedding theorem (Adams 1975). Inequality (3.17) is a consequence of equation (3.16) and the Schauder estimate for equation (2.3).

Proposition 2.1 combined with theorems 3.1 and 3.2 now yield important information about solution branches. For example, if a solution branch *Σ*_{n} is unbounded in , as described in proposition 2.1, then theorems 3.1 and 3.2 ensure that *Σ*_{n} has an unbounded projection onto . On the other hand, if only the second characterization of *Σ*_{n} holds as in proposition 2.1, *viz*. *Σ*_{n} meets the trivial line at another point , then *Σ*_{n} is bounded. In this case, we have the existence of solutions to equation (2.12) (for *λ*=1/2) in —either for all values *κ*_{n−1}<*κ*<*κ*_{n} or for *κ*_{n}<*κ*<*κ*_{n+1}. Either way, we have the existence of global solutions for arbitrarily large values of *κ* (possibly by taking *n* sufficiently large).

Similarly, we may fix *κ* sufficiently large leading to equation (2.11). Again, we have a general Rabinowitz-type result for the branches —now with *λ* as the parameter, cf. proposition 2.1. In this case, theorems 3.1 and 3.2 show that all solutions emanating from the trivial line (according to equation (2.11)) are contained inside a ‘cylinder’ , where *B*_{C} is a ball in , centred at the origin, of some radius *C*. In particular, these solution branches fulfil the second alternative of proposition 2.1, i.e. each branch ‘returns’ to the trivial line. This type of behaviour is similar to that observed in one-dimensional and two-dimensional Cahn–Hilliard models (e.g. Lifshitz & Rybnikov 1985; Triantafyllidis & Bardenhagen 1993; Kielhöfer 1997; Vainchtein *et al*. 1999; Maier-Paape & Miller 2002).

## 4. Numerical results

As discussed in remarks 3.1, the fruits of the previous two sections provide basic existence results and a specific characterization of the ‘location’ of global solution branches to our problem (2.10). We are now in a position to compute detailed solution structures—treating either *λ* or *κ* as the free continuation parameter. In particular, we are interested in solutions corresponding to local minima of the energy and large values of *κ*. The former terminology is characterized with precision shortly.

All results reported herewith were obtained by the following strategy: We first fix *λ*=1/2 (at the peak of the spinodal region, cf. (2.2) and figure 1), and then compute global solution branches via standard path following—with *κ* as the free parameter. Figure 2 depicts a sample global bifurcation diagram resulting from such computations. It includes parts of the continua bifurcating from the trivial solution at the bifurcation points (*κ*_{n}, 0), corresponding to the modes *u*_{n}=sin *nx* sin *y*, *n*=2, 3, …, 6, cf. (2.5) and (2.12). We call these *primary* branches. Here and in the sequel, we plot on the vertical axis. The black horizontal line on the bottom in figure 1 is the trivial solution *u*≡0. A few secondary bifurcation branches, also shown in figure 1, are not considered further in this work. Except for the primary branches emanating from the trivial line, the apparent intersections of the primary branches in figure 1 are not true, each being a consequence of the projection onto the plane.

We now continue the description of our solution strategy. We look for primary solution branches that are (apparently) unbounded. As mentioned in remarks 3.1, such branches contain solution pairs corresponding to arbitrarily large values of *κ*. We continue computing points on such branches until we find solutions corresponding to local minima of the potential energy, i.e. *stable* solutions. We now make the latter precise: first we rewrite (2.1) in terms of *u*=*w*−*λx* and *κ*=1/*ϵ*(4.1)For fixed values of *λ* and *κ*, we seek solutions, , at which the second variation (4.1) is positive definite, i.e.(4.2)for all , ** η**≠0, where .

We first discretize all derivatives in equation (2.3) via a finite-difference method. We then use a path-following algorithm to compute solution branches of the discretized system. A detailed description of the former can be found in Maier-Paape & Miller (2002). Here, we use proposition 2.1 to a great advantage: In computing each branch, , *n*=2, 3, …, it is enough to solve the problem on the subdomain —due to the inherent periodicity and oddness properties, cf. (2.9). Once a numerical solution for the subdomain problem (with Dirichlet boundary conditions) has been computed, the solution on the entire domain is easily generated by inversion–reflection (2.7) and 2*π*/*n*-periodic extension. Of course, the complete solution on *Ω* is required for the stability test (4.2), the latter of which is always carried out in . Moreover, the discretization employed in equation (4.2) must be consistent with that used to obtain the solution branch itself. For the latter, each branch is obtained using a uniformly rectangular grid on *Ω*_{n}. Results for grids of sizes 50×50, 75×75 and 100×100 were carried out for each primary solution to ensure reliability of the numerical approximation. We then check for stability in equation (4.2) using the results from the 100×100 grid on *Ω*_{n}. The aforementioned extension of the solution on the whole of *Ω* yields a (100*n*+*n*−1)×100 grid, i.e. the discretized vector ** η** in equation (4.2) has (100

*n*+

*n*−1)×100 entries. Accordingly, the local stability of a given computed solutions is tantamount to the positive-definiteness of the [(100

*n*+

*n*−1)×100] × [(100

*n*+

*n*−1)×100] Hessian matrix resulting from the finite-difference approximation of equation (4.2). That matrix is sparse and symmetric. This task is carried out efficiently by mostly tracking a small number of eigenvalues of the Hessian having values close to zero while proceeding along a solution branch.

More precisely, we employed the following strategy for determining stability: We first use a built-in Matlab function for sparse matrices, called *eigs*, to compute a few eigenvalues around 0. If we find a negative eigenvalue, we stop. Otherwise, we then use the Matlab program *eigifp*. For a symmetric matrix it returns the smallest eigenvalue, and it is possible to improve convergence and calculation time by using a suitable starting point. This is necessary in our case, and for this purpose, we use the smallest eigenvalue and its corresponding eigenvector computed by *eigs*. In any case, all solutions reported herewith as being stable had corresponding Hessians for which *eigifp* returned a positive eigenvalue.

To conclude the discussion of our solution strategy, once we find a locally stable solution on a *Σ*_{n} branch with a large-*κ* value, we fix that particular value of *κ* and then proceed with *λ*-continuation from the previously fixed value *λ*=1/2. As before, we can track the stability of solutions along such *λ*-branches.

The solution procedure just described delivered stable solutions, e.g. for *n*=7, 8, …, 12, with corresponding values of *κ* in the range 10 000–20 000. Figures 3–10 depict sample results from four of these calculations. Figures 3–6 show the global *λ*-branches for *n*=8, 9, 11 and 12, respectively. Observe that each is in keeping with remarks 3.1. Note the six points on each solution branch, denoted by small ‘o's’. In figures 7–10, we depict the associated solutions, all of which are stable. Figures 7–10 are contour plots of the associated solutions, which are constructed as follows: We use three shades for the level sets of the total shear strain *γ*_{x}≔*λ*+*u*_{x}, as encoded in figure 1. All regions of *Ω* on which the value of *γ*_{x}(*x*, *y*) lies within the ‘left well’ of are shown dark, and likewise those within the ‘right well’ are grey. Otherwise, for *γ*_{x}(*x*, *y*) within the so-called spinodal region, we choose the colour white, cf. figure 1. Note the relatively thin zones of white in figures 7–10. We emphasize that the solutions represented in each of the zones—dark, grey and white—in figures 7–10 are *not* homogeneous. Finally, we observe from equation (2.3) that (*κ*_{0}, 1−*λ*_{0}, −*u*_{0}) is a solution of equation (2.10) whenever (*κ*_{0}, *λ*_{0}, *u*_{0}) is a solution. Accordingly, each of the stable solution points denoted in figures 3–6 has such reflected solutions that are each stable and yield figures identical to those in figures 7–10, except that the grey and dark shades are reversed.

## 5. Concluding remarks

Our approach to finding locally stable equilibria corresponding to two-phase striped mixtures has several key ingredients: (i) first, we perform a global bifurcation analysis using the reciprocal of the capillarity coefficient, *κ*=1/*ϵ*, as the bifurcation parameter. In essence, we run through families of models (in the capillarity coefficient *ϵ*) seeking equilibria that correspond to small values of *ϵ*. Of course, the procedure of simply dividing the governing equation by *ϵ* alone is not enough to guarantee the existence of physically appropriate solutions—even an unbounded branch of solutions (according to the Rabinowitz alternative) need not contain large values of *κ*. (ii) Next, we establish *a priori* bounds. In this way, we know before computing that our global branches contain solutions corresponding to arbitrarily small capillarity. While this is suggestive of the existence of stable solutions, we need to compute next. (iii) Finally, we compute solution branches and explicitly check for stability of the discretized system.

We point out that we are solving a strongly elliptic problem (2.3) in fixed-point spaces (2.9). In particular, we tailor our finite-difference mesh to the symmetry inherent in the fixed-point space. Accordingly, we do not suffer from the problem of mesh-dependent solutions, as can be the case in the computation of minimizing sequences in sharp-interface theories (Luskin 1996) (*ϵ*=0 in our formulation).

Recall that each of the solution branches depicted in figures 3–6 were first determined by computing *κ*-branches, *Σ*_{n} (cf. proposition 2.1), until stable solutions were found. Next the *λ*-branches, , shown in figures 3–6 were computed (for fixed *κ*) from the ‘top down’, i.e. from the maximum norm solution shown down to near the *λ*-axis. It is natural to wonder why we did not simply pick large *κ*-values (say, by trial and error), use (2.11) and then directly compute the branches , in the hope of finding stable solutions. There are two essential reasons: First, that approach is tedious to say the least. Second, the numerical computations of the branches near the trivial solution in figures 3–6 are not reliable. Our method directly computes the upper branches, containing stable solutions, the computations of which are very reliable.

We only determine local stability of our solutions via the second variation, cf. (4.2). In particular, we do not know if our solutions correspond to global energy minima. Indeed, some of our solutions are probably not global minimizers. Note that *n*=8, 9 solutions in figures 3 and 4, respectively, each exist at *κ*=15 000. Similarly, *n*=11, 12 solutions in figures 5 and 6, respectively, each exist at *κ*=20 000. In each of these two instances, we have two distinct stable solutions for the same parameter values. It is highly unlikely that the two solutions correspond to the same energy value, i.e. at best only one of them would be the global energy minimizer.

Returning to figures 7–10, observe that the white zones are most visible near the top and bottom horizontal boundaries of the domain. Indeed, for each of these equilibria, the hard-loading boundaries are compatible with the large-shear grey phase but incompatible with the small-shear dark phase. In contrast, the vertical white zones between the grey and dark stripes are nearly imperceptible in the figures. Clearly, we have two different transition-boundary length-scales. This is strikingly similar to the two different length-scales commonly observed in real experiments—austenite and twinned-martensite boundaries on the one hand versus twin boundaries on the other (cf. Chu 1993).

In this work, we have exclusively focused on primary global solution branches. In particular, we have not pursued symmetry-breaking bifurcations. As hinted in §2, cf. (2.7), our (extended) problem has hidden O(2) symmetry, and each of our primary branches represent solutions with *D*_{n} symmetry for *n*=2, 3, …. It is conceivable, if not highly likely, that there are secondary bifurcations connecting commensurate primary branches, e.g. a branch connecting *n*=6 primary branch to *n*=12 branch. If so, then one would expect to find solutions characterizing the so-called ‘phase-tip splitting’ or ‘branching’ (cf. James *et al*. 1995) near the top and bottom horizontal boundaries of the primary domain—a phenomenon typically observed in experiments (Chu 1993). We plan to investigate this question in a future work.

Of course, our anti-plane shear problem considered here—although more complex than the Cahn–Hilliard model—is simpler than a bona fide two-dimensional elasticity model like those considered in Ball & James (1987), Ball & James 1992, and Bhattacharya (2003) (albeit in the absence interfacial energy). Our approach presented here holds great promise for the analysis of such problems in the presence of both interfacial energy and loading.

Although we do provide rigorous existence results in this paper, we point out that our numerical results presented herewith are not rigorous. However, there exist new methods that could potentially show that our computed equilibria and their stability are rigorous (Miller 2005). The basic idea is to embed our static problem within the framework of gradient dynamics. The time-dependent equation is written as an infinite dimensional system of ordinary differential equation using a complete orthogonal basis of a suitable Hilbert space. Conley index theory is then used to determine existence of a solution within a small neighbourhood of the solution obtained by a numerical approximation.

## Acknowledgments

We thank P. Rosakis and S. Maier-Pappe for useful suggestions. The work was supported in part by the National Science Foundation through grants DMS-0072514 and DMS-0406161, which is gratefully acknowledged.

## Footnotes

- Received October 3, 2006.
- Accepted December 14, 2006.

- © 2007 The Royal Society