Flexoelectricity, incommensurate phases and the Lifshitz point

The solutions for the minimizers of the energy density f (q, p)  =  Aq2 + Bq4 + p2 + gA,B + β(q′p−p′q) + |q′|2 + κ|p′|2 ?> describe the flexoelectric effect with a flexoelectric coupling coefficient β. The order parameters q and p can be visualized as strain and polarisation, respectively. The parameter κ denotes the ratio of intrinsic length scales for q and p. We show that the structural ground-states include 3 phases, namely the paraelastic state q  =  p  =  0, the ferroelastic state where polarization exists inside and near twin boundaries, and the incommensurate (modulated) phases with a very rich array of structural modulations ranging from nearly pure sine waves to kink arrays and ripple states. The phases coincide in the multicritical Lifshitz point. Linear flexoelectricity p∼q′ ?> is encountered only approximately inside the ferroelastic phase and near the phase boundary between the paraelastic phase and the incommensurate phase. The relationship between the polarisation and the strain gradient is highly non-linear in all other states. The spatial profiles and energy distributions are discussed in detail.

The solutions for the minimizers of the energy density f (q, p) = Aq 2 + Bq 4 + p 2 + g A,B + ( ) β κ − +| | + | | ′ ′ ′ ′ q p p q q p 2 2 describe the flexoelectric effect with a flexoelectric coupling coefficient β. The order parameters q and p can be visualized as strain and polarisation, respectively. The parameter κ denotes the ratio of intrinsic length scales for q and p.
We show that the structural ground-states include 3 phases, namely the paraelastic state q = p = 0, the ferroelastic state where polarization exists inside and near twin boundaries, and the incommensurate (modulated) phases with a very rich array of structural modulations ranging from nearly pure sine waves to kink arrays and ripple states. The phases coincide in the multicritical Lifshitz point.
Linear flexoelectricity p q ∼ ′ is encountered only approximately inside the ferroelastic phase and near the phase boundary between the paraelastic phase and the incommensurate phase. The relationship between the polarisation and the strain gradient is highly non-linear in all other states. The spatial profiles and energy distributions are discussed in detail.
Keywords: flexoelectricity, Lifshitz multicritical point, incommensurate phase, gradient coupling, order parameter coupling, polar twin boundary (Some figures may appear in colour only in the online journal) Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. linear-quadratic coupling for volume effects and incommensurations [27][28][29], and bilinear coupling in feldspars [30], and trilinear coupling [31] in PbZrO 3 require that the irreducible representations of order parameter products contain the identity representation. Such restrictions do not apply for the flexoelectric effect so that the unknown mix between various coupling schemes may be one reason why flexoelectricity in some materials far exceeds theoretical predictions [14,[32][33][34]. Polarization switching in incommensurate phases induced by flexo-or 'gradient-coupling' also exists in ferroelectric thin films [35] with little effect from mechanical boundary conditions [36]. Flexoelectricity represents hence a special case of coupling between order parameters and their gradients rather than between homogeneous order parameters [11,26,[37][38][39]. Such coupling between structural gradients and dipole moments not only generates local polarization but also renormalizes the gradient energy itself. The gradient energy is always positive in the uniform state but reduces by coupling. A transition between a uniform phase and a modulated (incommensurate) phase is reached when the gradient energy vanished. In addition, a multicritical Lifshitz point [38,[40][41][42][43][44][45][46][47][48] occurs when both the gradient energy and the uniform order parameter vanish simultaneously. It is the purpose of this paper to report that this mechanism is indeed an intrinsic property of the flexoelectric effect. We show that the gradient coupling thus generates a complex interplay of several structural phases where aperiodic patterns are separated from incommensurate phases by a multicritical Lifshitz point.
We consider a 'minimal' model with a one-dimensional energy density, which we construct initially from two order parameters Q and P. The motivation is to derive a minimum model with an anharmonic potential for Q and a purely harmonic energy density of P. The coupling is the gradient coupling between Q and P, their Ginzburg (bending) energies are allowed to be different. The low dimensionality of the energy density covers models in higher dimensions when the pattern consists of stripe patterns in the projection perpendicular to stripes in three dimensions or lines in two dimensions. We do not consider vortex patterns in this paper; they were discussed in [11] and [36]. We will show that the physics of the phase transitions is fully described by only one, scalar order parameter (the second order parameter relaxes into its energy minimum). The governing Euler-Lagrange equation contains only one order parameter. The phase diagram is the same as in [49] for a type II transition. For this reason we refer to the singular point in the phase diagram as 'Lifshitz point' rather than 'triple point' although we can derive only bounds for the trans itions. A detailed analysis of the criticality of the 'Lifshitz point' is beyond the purpose of this paper and will be taken up in a forthcoming contribution.

The total energy
We write the Gibbs free energy density of the uniform contribution as , : → R R represent Q and P, respectively, and f denotes the energy density. It suffices to consider positive β since for 0 β = the problem decouples, and for 0 β < the problem can be transformed to the case 0 β > by a sign change for either q or p. The ansatz for κ is not restricted to 1 κ = which would imply the singular situation where the two intrinsic length scales of the two order parameters are equal.
Since qp q p p q ( ) ′ = + ′ ′ is a null Lagrangian, the flexocoupling term in (1) can equivalently chosen to be q p 2 ′ or p q 2 − ′ . The corresponding system of Euler-Lagrange equations is Combining both equations yields the governing Euler-Lagrange equation: Bq q Bq q 1 2 12 6 . (3)

The minimization problem
We want to consider periodic and non-periodic pairs (q, p) with minimal energy. Since for periodic pairs E q p , 0, [ ] { } ∈ ±∞ , a refined notion of minimizer adapted from [50] will be applied. First, we define the energy growth rate E * of the total energy E:  Finally, if (q, p) is a locally minimal energy configuration and attains the minimal energy growth rate, we call it a minimal energy configuration.

If for some parameters A B
, , , β and κ, the minimal energy rate is non-negative, then we will consider the minimizers (in the classical sense) of E with boundary conditions as the solution to our minimization problem and call them minimizers. These are the typical kinks and breathers connecting the zeroes of G(Q, P), and they always have energy rate 0. If the minimal energy rate is negative, we will consider a minimal energy configuration as a minimizer. For the energy at hand, these are always periodic and, in particular, globally bounded.

Minimizers and the Euler-Lagrange equations
The Euler-Lagrange equation (3) belongs to a class of equations that has a very rich set of bounded solutions. A model equation for this class is which is satisfied by stationary solutions of the extended Fisher-Kolmogorov equation for 0 λ < and stationary solutions of the Swift-Hohenberg equation for 0 λ > . For 8 λ > − , the solution space of (5) consists of an infinite number of kinks and pulses with arbitrarily many 'bumps' as well as periodic solutions that can be assembled out of kinkand pulse-like bits [51].
Regarding the total energy or the energy growth rate, these solutions are at most metastable since being a minimizer implies higher symmetry constraints than solving the Euler-Lagrange equations. For a discussion of the stability of such solutions see e.g. [52,53].
For a physics perspective, the two approaches represent different scenarios. In case of the Euler-Lagrange equation and the equivalent pattern formation approaches one describes the transformation of an, often uniform, material into a structured pattern. The transformation progresses from a nucleation centre whereby the mechanism is determined by the behaviour of the propagating front [54,55]. In contrast, the energy minimizer determines the pattern which has the lowest energy. In many situations the structures are sufficiently equilibrated that only the pattern with the lowest energy is expected to be observable. This is the approach taken in this paper.

Results
First, we discuss the case B > 0. Afterwards, we treat the case B = 0 as an approximation of the regime B  A.

The case B > 0
For B > 0, the phase diagram spanned by A ∈ R and 0 β > is visualized in figure 1 and consists of three phases: • a paraelastic phase were the minimizers are zero, • a ferroelastic phase where the minimizers are kinks and breathers, • an incommensurate phase where the minimizers are periodic.
The incommensurate phase consists of the parameters for which the minimal energy growth rate is negative. In the para-and ferroelastic phase, the minimal energy growth rate is non-negative and, thus, the minimizers (q, p) are classical minimizers of E[q, p] that satisfy the boundary conditions (4). If A 0 ⩾ , these minimizers connect the point (Q, P) = (0, 0) with itself, wherefore they are zero themselves, and the parameters belong to the paraelastic phase. If A < 0, the minimizers connect the points (±Q 0 , 0), and the parameters belong to the ferroelastic phase.
The boundary between the paraelastic phase and the incommensurate phase is independent of B and given by This relation will be derived in the sections 4.1.1 and 4.1.2 below.
The boundary between the ferroelastic phase and the incommensurate phase is implicitly given by a critical , ,κ the minimal energy growth rate changes sign at β * . We are unaware of a more explicit description. However, a lower bound ( 1 β > * ) and an upper bound (13) are derived below. The numerically obtained boundaries for three values of κ are visualized in figure 2.
The Lifshitz point at which the three phases touch is independent of B and κ and located at A 0, 1 β = = . The phase lines in figure 1 approach each other tangentially, as predicted for a Lifshitz point of type II in [49].

Lower bounds for the incommensurate transition.
To obtain lower bounds in β for the boundaries to the  incommensurate phase we estimate the energy growth rate E * from below. To this end, we use that that pq p q q p ( ) ′ = + ′ ′ is a null Lagrangian and that the periodic minimizers are bounded absolutely by some constant M > 0, wherefore where we have dropped non-negative terms, used the null Lagrangian, and applied Cauchy's inequality, i.e.
with 1 ε = . Thus, the minimizer (q, p) is a classic minimizer of E[q, p] that connects the points (±Q 0 , 0). If 1 β > , the argument can be amended with the additional constraint Then, using similar arguments as above yields has been chosen in the last step. In the paraelastic phase, i.e. A 0 ⩾ , the bound (7) corresponds to the boundary (6) given above. In the ferroelastic phase, i.e. A < 0, (7) cannot be satisfied wherefore we only have the lower bound 1 β > * . Note that both estimates above are independent of B since the term Bq 4 was simply dropped.

Upper bounds for the incommensurate transition.
To obtain upper bounds in β for the boundaries to the incommensurate phase, we calculate the energy growth rate for trial functions and check for which parameters it is negative. Let a, b, k > 0 and consider the ansatz Then, Optimizing the parameters a, b, k under the constraints yields the parameters and that the minimal energy growth rate for the ansatz is The details of the calculation are given in the appendix.
Since g A, B = 0 for A 0 ⩾ , it follows together with the discussion from section 4.1.1 that (6) is indeed the boundary between the paraelastic and the incommensurate phase. Numerical simulations show that the periodic minimizers near the boundary are very similar to pairs of sine and cosine (figure 3, Case (i)). While the similarity grows stronger with decreasing distance to the boundary, the amplitudes of the minimizer vanish as indicated by (11). Note that (11)    Together with the lower bound from the discussion above, this implies: In figure 2, the left coordinate axis and the blue line mark the lower and upper bound, respectively. The transition lines from numerical simulations show that the upper bound is good for small κ, but rather rough for large κ. Hence, the ansatz (8) resembles the minimizers along the phase boundaries only well for small κ.
Note that at the Lifshitz point the modulation wavevector k is always zero as required for a type II behavior. Figure 3 shows plots of the minimizers for representative parameters and figure 4 shows the corresponding relations q p ( ) ′ and p q ( ) − ′ . The four cases differ in the values of A and β, while B = 1 and 5 κ = are fixed:

Qualitative behaviour of minimizers.
The parameters belong to the incommensurate phase and are near the transition to the paraelastic phase. As suggested by the discussion above, the minimizer  The parameters belong to the incommensurate phase and β is large. Over longer distances, q is almost constant and p behaves almost linearly. These distances are connected by kinks in q and smooth extrema in p. While the relation q p ( ) ′ is almost linear, the relation p q ( ) ′ is again highly non-linear.
In both the ferroelastic and the incommensurate phase the amplitudes of the minimizers converge to zero as the Lifshitz point is approached. All plots shown are the result of numerical simulations whose core routine was the minimization of the total energy (1) on finite intervals. The periodic minimizers have been obtained by using periodic boundary conditions and an external minimization procedure to find the minimal energy growth rate.

The case B = 0
In the limit case B = 0, we restrict ourselves to A > 0 to ensure that G Q P , 0 ( )⩾ . Thus, there are only a paraelastic and an incommensurate phase.
Because of the scaling the minimization problem degenerates since the minimal energy growth rate is either 0 in the paraelastic phase or −∞ in the incommensurate phase. We will show below that the boundary between the two phases is again given by (6) and that in the incommensurate phase the space of bounded solutions to the Euler-Lagrange equations is spanned by pairs of sines and cosines. The estimates in section 4.1.1 only required B to be nonnegative. Therefore, we already have the lower bound (7) for the boundary between the two phases.
For the upper bound we use again the ansatz (8) and obtain with b and k as in (11) , then the growth rate becomes arbitrarily small with growing a, which shows the upper bound.
Setting B = 0 in (3) yields the governing Euler-Lagrange equation: This equation has bounded non-vanishing solutions if and only if its characteristic polynomial has (purely) imaginary roots. Since (16) is a parabolic equation with respect to k 2 and since the lowest-order coefficient is positive, the requirement is satisfied if the coefficient of k 2 and the discriminant D are non-negative. In the incommensurate phase, this is true since 1 where k 1 and k 2 are the two imaginary roots of (16). Substituting these solutions into (2) yields the equations p x p x q x a k k x cos .
Hence, for the two basis solutions the relations q p ( ) ′ and p q ( ) ′ are linear.

Discussion
Flexoelectricity depends on three parameters: the entropy term of the driving order parameter A, the flexoelectric coefficient β and the ratio of the length scales κ. The 'classic' flexoelectric effect between two order parameters with parabolic selfenergies exists for A ≫ B with the solution B = 0 for systems without phase transitions. The ground state is then paraelastic and the linear flexoelectric relation p q = ′ (and q p = ′) is always fulfilled. This situation is encountered when a sample x f is bent macroscopically and induces surface charges [25]. For B > 0 a ferroelastic phase transition occurs for A = 0 and 1 β < . The Lifshitz point is situated at the end of this interval at A = 0 and 1 β = . The Lifshitz point is independent of κ and agrees with the requirements of a type II behavior in the nomenclature of [49].
A more detailed analysis of the asymptotic behavior near the Lifshitz point (or the triple point in a type I behavior [56]) will be presented in a forthcoming paper. Very careful experimental observations of the soft optic and acoustic branches near a type II Lifshitz point have been investigated in [57].
The ferroelastic phase transition at A = 0 and 1 β < leads to twinning under the appropriate boundary conditions. Polarity emerges from the twin boundaries as the desired key feature of domain boundary engineering [1] (figures 1 and 2). The order parameter diagrams are approximately semicircles in close analogy with the solutions for bilinear coupling [27]. Experimentally it will be very hard to distinguish between these two coupling schemes based on the polar deformation patterns inside the twin walls. Computer simulations of a simple model confirm this behaviour [36].
The third structural state is the incommensurate phase with 1 β > and is situated beyond the Lifshitz point. We denote this state 'incommensurate' (and not just modulated) because no lattice interaction exists so that the periodicity is unrelated to the lattice repetition length and phason excitations possess zero energy. This structural state is characteristic for the flexoelectric coupling and does not exist in any of the other coupling schemes. The proximity of incommensurate phases in a phase diagram is hence always an indication that flexo-polarity exists in these systems. A compilation of modulated crystal structures, where flexoelectricity is expected to play a key role for the structural modulation, was recently derived [58].
We now discuss a structural evolution with fixed ratio of length scales κ and constant flexoelectric coefficient 1 β > . Lowering A from the paraelastic state at A ≫ B, a transition to an incommensurate phase occurs at . The structural state near the transition point is characterized by harmonic waves in q and p and the flexoelectric relation p q = ′. This represents a 'classic' incommensurate phase with sin-wave modulations, which evolves into a soliton lattice at lower values of A. The spatial variation of q becomes a sequence of kinks while p shows almost triangular profiles (figure 1). The overshoots in the q profile near the centre of kinks are systematic; they are a consequence of the reduction of the energy density in this regime (figure 5). Such squared-up periodic patterns can change to aperiodic patterns if β is reduced. Near the transition points, which now depend strongly on κ, a multitude of metastable states form characterized by ripples, which are decaying oscillations pinned at the centre of kinks [11,37]. This transition regime is particularly rich in metastabilities and transient patterns with slow kinetics.
The role of the two length scales is different for the two transitions between the incommensurate phase and either the paraelastic phase or the ferroelastic phase. The transition between the incommensurate and paraelastic phase depends only on 1  . Reducing κ enhances the tendency for modulations ( figure 2).
The flexoelectric correlation p q ( ) ′ is approximately linear in the paraelastic phase and close to the transition between the paraelastic and the incommensurate phase. No linearity exists in other regimes where p increases steeply near q 0 = ′ and becomes weakly dependent on q′ for larger values of q′ ( figure 4). These S-shaped profiles reflect the large plateau in q between kinks where q′ is small. The polarization changes strongly in the same part of the pattern. This S-shape becomes well developed for 1 β and A = 0. In this case, we find no overshoots in the q profiles while p changes almost linearly over the entire period of the modulation. The order parameter diagram becomes almost identical to that of the bi-quadratic coupling in the strong coupling limit (as is the circle in case of the proximity of the paraelastic-incommensurate transition). This demonstrates that any quantitative analysis of polar patterns for the distinction between different coupling mech anisms is virtually impossible in these cases. The characteristic feature for flexoelectric coupling is the overshoot in q, which leads to pumpkin-shaped diagrams in figure 3. Since all variables are positive and 1 ⩾ β , rearranging the terms gives the unique solution b a 1 . That this critical value corresponds indeed to a minimum of E * , follows from its uniqueness, the fact that C is non-negative for b 0 → and b → ∞, and that finally C will be non-positive as can be seen from a comparison of (A.2) and (A.4). To find the representations of b and k given in (11) as in (11). Note that the right hand side is non-negative because of the assumption (10).