Coexistence of stable branched patterns in anisotropic inhomogeneous systems

Pattern formation is one of the most fascinating and intriguing phenomena in nature [1, 2]. It takes place in a wide variety of physical, chemical and biological systems and on disparate spatial and temporal scales, for example, convection phenomena in geoscience [2, 3] or in liquid crystals [4–6], environmental patterns [7–9], or patterns in chemical reactions [10, 11] and bacterial colonies [12]. In some circumstances pattern formation is undesired, for instance, the formation of spiral waves leading to cardiac arrhythmias in the heart muscle [13]. In other contexts pattern formation is even essential for the functioning of a system, e.g. in embryo development [14] or when designing surface wrinkling patterns to fabricate nanometer-scale structures [15, 16]. The mechanisms leading to the same type of pattern in different systems are obviously very diverse. Nevertheless, patterns occuring in systems of the same symmetry share common qualitative properties that can be described by universal amplitude equations for the envelope of periodic patterns [2, 5, 6, 17].

Generating, modifying or eliminating patterns hence either requires a profound understanding of the pattern formation mechanism in each specific system, or complementary, of the universal properties of a class of patterns and their response to symmetry breakings. Common pattern interventions are feed back control [11] and symmetry breaking via spatial, temporal or spatio-temporal modulations [19–30]. To mention two interesting scenarios, spatial forcing near resonance (between the forcing and the natural wavelength) can lead to so-called incommensurate patterns [19, 20], while symmetry breaking via long-wave spatial modulations can render stationary patterns time dependent [31–34]. The response of patterns in quasi one-dimensional (1D) systems is now well established, however truely 2D scenarios, like the interplay of an anisotropy and a modulation in different directions, are yet fully unexplored.

In this paper we identify and analyze a new class of quasi-2D anisotropic inhomogeneous pattern forming systems. The wave vector of the patterns lies close to ${\hat{{\bf{q}}}}_{0}=({\hat{q}}_{0},0)$ along the preferred x-direction (anisotropy), as for the two experimental systems sketched in figure 1: (a) a spatially modulated version of the wrinkle forming system [16, 35] and (b) a modulated version of electroconvection (EHC) in nematic liquid crystals [6]. By a modulation we break in such systems the translational symmetry along the perpendicular y-direction, causing a variation of the pattern's natural wave number ${\hat{q}}_{0}.$ This can be accomplished by varying the elasticity in the wrinkle forming system or the height of the electroconvection cell, respectively. In this class of anisotropic systems we find straight stripes, see figure 1(c), that are stable for wave numbers in a finite range around ${\hat{{\rm{q}}}}_{0},$ similar as in homogeneous anisotropic systems. In addition, however, we find a whole family of stable branched patterns as shown in figure 1(d) and (e). Surprisingly, they have—at identical parameters—different characteristic wave numbers and include different numbers of branching points. Moreover, the branched patterns coexist with the straight stripes in a wide parameter range. This behavior is a non-trivial generalization of the wave number bands (Eckhaus bands) for homogeneous systems [37–43]—a well established concept and experimentally verified e.g. in EHC [39] and axisymmetric Taylor vortex flow [41, 42]—to inhomogeneous systems and multiple patterns. In the following we present the universal amplitude equation of this new symmetry class of patterns and analyze its solutions.

Zoom In Zoom Out Reset image size

A generic model for the formation of stationary periodic patterns in anisotropic—but homogeneous—2D systems, described by a field $u(x,y,t),$ has been proposed in [43]. One interpretation of the field $u(x,y,t)$ is, that it describes the (small) lateral displacement of a thin elastic plate extended in the $x-y$ plane, loaded along the x and y direction and supported by an elastic medium, similar as in wrinkling systems (whereby the in-plane elastic deformations are neglected) [43]. Here we generalize this model to an inhomogeneous situation by modulating the pattern's preferred natural wave number ${\hat{{q}}}_{0}$ along the direction perpendicular to the anisotropy (the x-direction),

$$\begin{eqnarray}&&{q}_{0}(y)={\hat{q}}_{0}+M\mathrm{cos}({k}_{m}y),\end{eqnarray} \tag{ 1 }$$

with an amplitude M and a wave number k_m considerably smaller than ${\hat{q}}_{0}.$ Then the dynamics of the patterns, described by the real field $u(x,y,t)$ is governed by

$$\begin{eqnarray}{\partial }_{t}u & = & \left[\varepsilon -{\left({q}_{0}^{2}(y)+{\rm{\Delta }}\right)}^{2}\right]u-W{\partial }_{x}^{2}{\partial }_{y}^{2}u-c{\partial }_{y}^{4}u-{u}^{3}\\ & & -2({\partial }_{y}u){\partial }_{y}\left[{q}_{0}^{2}(y)\right]-u{\partial }_{y}^{2}\left[{q}_{0}^{2}(y)\right].\end{eqnarray} \tag{ 2 }$$

The first line corresponds to the original model in [43] and the second line includes contributions due to the modulation ${q}_{0}(y).$ Equation (2) is a representative of the here-identified symmetry class. It can be directly linked to the elastic wrinkle-forming system via an appropriate rescaling [43]: ${q}_{0}^{4}=\frac{\kappa }{{\lambda }_{1}}$ relates the critical wave number to the bending stiffness λ₁ of the hard layer and the elastic modulus κ of the substrate. The control parameter $\varepsilon =\left(1-\displaystyle \frac{{\mu }_{1,c}^{2}}{{\mu }_{1}^{2}}\right)\;{q}_{0}^{4}$ is related to the critical compression ${\mu }_{1,c}=2\sqrt{\kappa {\lambda }_{1}}.$ Also note that, similar to the homogeneous version [43] and the wrinkle system, equation (2) can be derived from a functional as described in the appendix. For the following, we have chosen W = 1, c = 0.5 and ${\hat{q}}_{0}=1$ that favor straight stripe formation parallel to the y-direction for the homogeneous case.

Close to the onset of supercritically bifurcating patterns their generic (system-independent) properties may be described by a nonlinear dynamical equation for the complex envelope $A(x,y,t),$ considering $u(x,y,t)=A{{\rm{e}}}^{{\rm{i}}{\hat{q}}_{0}x}+{A}^{*}{{\rm{e}}}^{-{\rm{i}}{\hat{q}}_{0}x}$ with the scaling $u\propto {\varepsilon }^{1/2}$ (the star denotes the complex conjugate). This reduction method, the so-called multiple scale analysis, is well established for supercritical bifurcations in 2D homogeneous isotropic [2, 18] and anisotropic systems [5, 43] or 1D inhomogeneous ones [32]. Following very closely the appendix in [2] with the same intermediate scaling for time and space one obtains with $M\propto {\varepsilon }^{1/2}$ and ${k}_{m}\propto {\varepsilon }^{1/2}$ via the multiple scale analysis from the system described by equation (2) the generic amplitude equation for the here-identified universality class of patterns in anisotropic systems with a spatially varying natural wave number ${q}_{0}(y):$

$$\begin{eqnarray}{\partial }_{t}A & = & \left[\varepsilon +{\xi }_{x}^{2}{\partial }_{x}^{2}+{\xi }_{y}^{2}{\partial }_{y}^{2}\right]A-g{| A| }^{2}A\;\\ & & -{\rm{i}}8{\hat{q}}_{0}^{2}\;M\mathrm{cos}({k}_{m}y){\partial }_{x}A-4{\hat{q}}_{0}^{2}\;{M}^{2}{\mathrm{cos}}^{2}({k}_{m}y)A.\end{eqnarray} \tag{ 3 }$$

Upon derivation from equation (2) one obtains the relations ${\xi }_{x}^{2}=4{\hat{q}}_{0}^{2},{\xi }_{y}^{2}=W{\hat{q}}_{0}^{2}$ and g = 3. The generic phenomenon of stripes' branching, as detailed below, is induced in equation (3) by the term $\propto M{\partial }_{x}A.$ The term $\propto {M}^{2}A$ mainly leads to quantitative modifications. It should be noted that the first line of equation (3) is exactly the equation for the amplitude of stripe patterns in unmodulated EHC [5, 43], the second anisotropic system suggested for experimental investigations of phenomena described here. For all specific pattern forming systems in the considered universality class, equation (3) can be derived via perturbation techniques from the basic equations (as for example from equation (2)) and it has always the same form—only the coefficients will depend on the specific system. The presented results are obtained via the universal equation (3).

The modulation of the natural wave number, see equation (1), alters the bifurcation from the basic state (u = A = 0) towards stationary periodic solutions of the form $A(x,y)={{\rm{e}}}^{{\rm{i}}{Qx}}H(y).$ The stability of the latter and of A = 0 with respect to small perturbations $v(x,y,t)={v}_{1}(y){{\rm{e}}}^{\sigma t+{\rm{i}}{Kx}}+{v}_{2}^{*}(y){{\rm{e}}}^{{\sigma }^{*}t-{\rm{i}}{Kx}}$ can be determined via the ansatz $A(x,y,t)={{\rm{e}}}^{{\rm{i}}{Qx}}\left[H(y)+v(x,y,t)\right],$ followed by a linearization of equation (3) with respect to small v. An expansion of H(y) and ${v}_{\mathrm{1,2}}(y)$ in Fourier modes results in an eigenvalue problem for σ, whereby solutions are stable when all eigenvalues have negative real parts. Note that in $v(x,y,t)$ we consider only perturbations along the x-direction, as perturbations parallel to the stripes (along y) have been numerically found to have no effect.

For the basic state A = 0 the stability condition, ${\mathfrak{R}}\left[{\sigma }_{\mathrm{max}}(\varepsilon ,Q)\right]=0,$ determines the neutral curve ${\varepsilon }_{{\rm{N}}}(Q),$ above which stationary periodic solutions, $A={{\rm{e}}}^{{\rm{i}}{Qx}}H(y),$ of finite amplitude do exist. The neutral curves for the unmodulated case M = 0 and for a small modulation M = 0.03 are plotted in figure 2(a) with dashed and solid lines, respectively. Both curves have their minimum at Q = 0 and are symmetric. However, when increasing M, the value of the minimum ${\varepsilon }_{{\rm{N}}}(0)$ shifts towards positive values due to the term $\propto {M}^{2}A$ in equation (3) that becomes larger ($\varepsilon \to \left[\varepsilon -4{\hat{q}}_{0}^{2}\;{M}^{2}{\mathrm{cos}}^{2}({k}_{m}y)\right]$). By increasing M further beyond a critical modulation (here M_c = 0.0459 for ${\hat{q}}_{0}=1$ and k_m = 0.1) a very important phenomenon takes place: The neutral curve develops two minima at finite wave numbers $Q=\pm {Q}_{\mathrm{min}},$ as shown in figure 2(b) by the solid line. The appearence of these two minima is caused by high modulation amplitudes and mainly by the term $\propto M{\partial }_{x}A$ in equation (3). The dependence of ${Q}_{\mathrm{min}}\propto \sqrt{M-{M}_{{\rm{c}}}}$ on the modulation in the vicinity of M_c resembles a pitchfork bifurcation, see figure 2(c). The two new terms in the amplitude equation, caused by the modulation, hence widen the neutral curve at small modulations and change its shape from a single- to double-well shape via a pitchfork-like bifurcation at higher modulations. The fact that the modulation along the y-direction causes two minima of the neutral curve ${\varepsilon }_{{\rm{N}}}(Q),$ with respect to Q along the x direction, is crucial for the emergence of branched patterns.

Zoom In Zoom Out Reset image size

Periodic solutions in the homogeneous case, M = 0, have a constant amplitude, $H(y)=\sqrt{(\varepsilon -{\xi }_{x}^{2}{Q}^{2})/g},$ and they are linearly stable only above the dotted line in figure 2(a), the so-called Eckhaus-stability boundary ${\varepsilon }_{{\rm{E}}}(Q).$ Note that for 2D anisotropic systems, the Q-range of stable solutions is symmetric with respect to Q = 0 [43], in contrast to isotropic systems, in which the zig-zag instability occurs for $Q\lt 0.$ For small modulations $M\lt {M}_{{\rm{c}}},$ H(y) is 2π/k_m-periodic along the y-direction and corresponds to modulated stripes as shown in figures 1(a) and 3(a), (b). In general, the equation for H(y) may have further solutions as well. However, simulations did not show them, so they are either unstable or have higher energies. The Eckhaus stability range of the modulated stripes is qualitatively unchanged compared to the homogeneous limit, as shown by the red-colored region in figure 2(a): it still touches ${\varepsilon }_{{\rm{N}}}(Q)$ at Q = 0. However, quantitatively, the modulation tends to narrow the width of the Eckhaus-stability range and thus reduces the range of stable stripes.

Zoom In Zoom Out Reset image size

Concomitant with increasing the modulation beyond the threshold M_c a striking change in the stability scenario occurs: (i) as shown in figure 2(b), a gap opens between the neutral curve (solid line) and the stability range of amplitude modulated straight stripes (red-colored area). Within the gap, the stripe solutions $A(x,y)={{\rm{e}}}^{{\rm{i}}{Qx}}H(y)$ are unstable and thus cannot be observed anymore. The width of this gap increases nonlinearly with M − M_c, as shown in figure 2(d). (ii) Stable branched stripes, like those shown in figures 1(d), (e) and figure 4, emerge within a large area above the neutral curve. (iii) Within the red-colored area in figure 2(b), the unbranched modulated stripes coexist with branched stripes at identical parameters sets.

Zoom In Zoom Out Reset image size

To characterize the coexistence (the first central result) further, we use the fact that the dynamical equations (2) and (3) can be derived from functionals, which for stationary patterns take the simple form (see the appendix):

$$\begin{eqnarray}&&{\mathcal{F}}=-\displaystyle \frac{1}{4}\int {\rm{d}}x{\rm{d}}y\;{u}^{4}=-\displaystyle \frac{3}{2}\int {\rm{d}}x{\rm{d}}y{| A| }^{4}.\end{eqnarray} \tag{ 4 }$$

Using the energy density $f={\mathcal{F}}/{L}_{x}{L}_{y}$ (${L}_{x}{L}_{y}$ is the domain area) one can hence determine which of the stable patterns is energetically preferred at a given parameter set (see also p. 868 of [2]). We have calculated the energy densities of unbranched patterns f_s and of branched patterns f_d within the red-colored region in figure 2(b). Their difference $\delta f={f}_{s}-{f}_{d}$ is plotted in figure 3(c) as a function of and for three different wave numbers Q. For small , branched patterns have a lower energy and thus are energetically preferred, whereas for considerably larger , unbranched patterns are in turn preferred. The crossover can be understood as the modulation becomes less important at large values of the control parameter .

The second central result is that for each parameter set out of a large range beyond ${\varepsilon }_{{\rm{N}}}(Q)$ we find always a whole family of stable branched patterns, characterized by different wavelengths and distinct numbers of branching points—and obtained by varying the initial conditions. For example, the two patterns in parts (d) and (e) in figure 1 contain five and seven pairs of branching points at identical parameters. They are both stable, however, the value of the functional of every branched pattern is, in general, different. This pattern coexistence in inhomogeneous systems is a surprising generalization of the well established Eckhaus wave number band of stable stripe patterns in homogeneous systems [37–42]. There, stable stripes of different wave numbers coexist at an identical parameter set, as verified in EHC [39], buckling [40] and axisymmetric Taylor vortex flow [41, 42]. Here, a band of branched patterns with different wavelengths and numbers of branching points coexists. In both cases, the origin of this coexistence lies in the fact that transitions between two stable periodic patterns require rather high excitations [38] worth to be investigated further.

Another interesting observation is due to the branching points (defects) having opposite topological charges in each half period π/k_m along the y-direction hence they repel each other. As shown in figure 4, this repulsion can lead to different orderings of the defects with increasing defect density, for instance, the zig-zag ordering shown in figure 4(b). It can even lead to more complex cascade-like ordering, as shown in figure 4(c). However, how the detailed ordering of branching points can be controlled by the magnitude, the wavelength and the anharmonicity of the modulation is an interesting and fundamental question that needs to be adressed in the future.

The phenomena identified here are robust and insensitive to the special shape of the imposed modulation of the natural pattern wave number. Studying the generic model, equation (2), and imposing instead of a harmonic modulation either an anharmonic wave-number modulation or a step-wise one, we obtain qualitatively the same results [35]. Consequently, instead of harmonic variations, which may be challenging to implement in experimental systems, a step-like modulation can be used, which is easy to achieve in the wrinkle system by glueing two elastomer substrates with different elastic properties together [35, 44].

In conclusion, we identified a new symmetry class of pattern forming systems that are anisotropic with a perpendicular modulation. By deriving and analyzing the universal amplitude equation, we showed that this class displays interesting scenarios—the emergence of families of stable branched patterns and their coexistence with unbranched patterns—that we suggest to verify experimentally in two complementary anisotropic systems: first, in wrinkle forming systems [16], where the elasticity of the substrate (as sketched in figure 1(a)) and/or the thickness of the hard layer on top can be varied perpendicularly to the pre-stretch (anisotropy) direction, as recently realized experimentally [35, 44]. And second, in dissipative electroconvection in nematic liquid crystals [4–6], where the layer height (as sketched in figure 1(b)) or the driving frequency are modulated perpendicularly to the nematic alignment. For both systems, a thorough theoretical analysis of the modulated basic equations is feasible, including the derivation of the universal equation of the patterns envelope, i.e. equation (3). There is a vast literature on wrinkles recently (see e. g. [15, 16, 35, 45–47]), including control strategies for wrinkle formation, but the question of pattern coexistence seems to have not even been raised yet. Hence studies along the lines proposed here, including quantitative comparisons between theory and experiments, will prove powerful for future control and design strategies, both of unbranched and branched (wrinkle) patterns. Finally, note that the branching of stripe patterns on the skin of fishes [48], an important example for patterns in living systems, is probably also induced by inhomogeneities, similar as discussed here, and which are in addition slowly changing with time.

We thank R Aichele, A Fery, W Pesch (Bayreuth), as well as V Delev (Ufa, Russia) for interesting discussions.

The dynamical equation (2) for the field $u(x,y,t)$ can be derived from the functional

$$\begin{eqnarray}{\mathcal{F}}[u(x,y,t)] & = & \displaystyle \int \tilde{f}(u){\rm{d}}x\;{\rm{d}}y=\displaystyle \frac{1}{2}\displaystyle \int \left[-\epsilon {u}^{2}+\displaystyle \frac{1}{2}{u}^{4}+W{\left({\partial }_{x}{\partial }_{y}u\right)}^{2}\right.\\ & & \left.+c{\left({\partial }_{y}^{2}u\right)}^{2}+\left({q}_{0}^{4}(y){u}^{2}+2{q}_{0}^{2}(y)u{\rm{\Delta }}u+{({\rm{\Delta }}u)}^{2}\right)\right]{\rm{d}}x\;{\rm{d}}y\end{eqnarray} \tag{ A.1 }$$

via the functional derivative

$$\begin{eqnarray}&&{\partial }_{t}u(x,y,t)=-\displaystyle \frac{\delta {\mathcal{F}}}{\delta u}=-\left(\displaystyle \frac{\partial \tilde{f}}{\partial u}+{\partial }_{x}^{2}\displaystyle \frac{\partial \tilde{f}}{\partial ({\partial }_{x}^{2}u)}+{\partial }_{x}{\partial }_{y}\displaystyle \frac{\partial \tilde{f}}{\partial ({\partial }_{x}{\partial }_{y}u)}+{\partial }_{y}^{2}\displaystyle \frac{\partial \tilde{f}}{\partial ({\partial }_{y}^{2}u)}\;\right).\end{eqnarray} \tag{ A.2 }$$

The functional in equation (A.1) can be simplified by assuming periodic boundary conditions and by using the expression ${\partial }_{t}u$ from the dynamical equation. One gets

$$\begin{eqnarray}&&{\mathcal{F}}[u(x,y,t)]=\displaystyle \frac{1}{2}\int \left[\displaystyle \frac{1}{2}{u}^{4}-u({\partial }_{t}u)-{u}^{4}-2u({\partial }_{y}u){\partial }_{y}({q}_{0}^{2}(y))-{u}^{2}{\partial }_{y}^{2}\left[{q}_{0}^{2}(y)\right]\right]{\rm{d}}x\;{\rm{d}}y,\end{eqnarray} \tag{ A.3 }$$

where the $-{u}^{4}$-term compensates for the cubic term occurring in ${\partial }_{t}u$ and similarly the last two terms. These last two terms again vanish due to the periodic boundary condition (here in y-direction):

$$\begin{eqnarray}&&\displaystyle \int \left[-2u({\partial }_{y}u){\partial }_{y}({q}_{0}^{2}(y))-{u}^{2}{\partial }_{y}^{2}\left[{q}_{0}^{2}(y)\right]\right]{\rm{d}}x\;{\rm{d}}y\quad =-\displaystyle \int {\partial }_{y}\left[{u}^{2}{\partial }_{y}\left({q}_{0}^{2}(y)\right)\right]{\rm{d}}x\;{\rm{d}}y=0.\end{eqnarray} \tag{ A.4 }$$

We arrive at the simple expression

$$\begin{eqnarray}&&{\mathcal{F}}[u(x,y,t)]=-\displaystyle \frac{1}{2}\int \left[\displaystyle \frac{1}{2}{u}^{4}+u({\partial }_{t}u)\right]{\rm{d}}x\;{\rm{d}}y\end{eqnarray} \tag{ A.5 }$$

and if we are—as in our work—interested only in stationary patterns, we simply have

$$\begin{eqnarray}&&{\mathcal{F}}[u(x,y,t)]=-\displaystyle \frac{1}{4}\int {u}^{4}{\rm{d}}x\;{\rm{d}}y.\end{eqnarray} \tag{ A.6 }$$

Also the amplitude equation (3) can be derived from a functional, which can be simplified for stationary solutions in the same manner, only the prefactor changes:

$$\begin{eqnarray}&&{\mathcal{F}}[u(x,y,t)]=-\displaystyle \frac{3}{2}\int | A{| }^{4}{\rm{d}}x\;{\rm{d}}y.\end{eqnarray} \tag{ A.7 }$$