Periodic boundary conditions for the simulation of 3D domain patterns in tetragonal ferroelectric material

Periodic domain patterns in tetragonal ferroelectrics are explored using a phase field model calibrated for barium titanate. In this context, we discuss the standard periodic boundary condition and introduce the concept of reverse periodic boundary conditions. Both concepts allow the assembly of cubic cells in accordance with mechanical and electrical conditions. However, application of the reverse periodic boundary condition is due to an increased size of the RVE and enforces more complex structures compared to the standard condition. This may be of particular interest for other multiphysics simulations. Additionally, we formulate mechanical side conditions with minimal spherical (hydrostatic) stress, or conditions with controlled average strain. It is found that in sufficiently small periodic cells, only a uniform single domain, or the simplest stripe domains constitute equilibrium states. However, once the periodic cells are of order 20 domain wall widths in size, more complex, 3-dimensional patterns emerge. Some of these patterns are known from prior studies, but we also identify other domain patterns with long, ribbon-like domains threaded through them and some vortex-like structures.

According to Maxwell's equations, domains usually arrange with electrically compatible 'head-to-tail' polarization, which minimizes energy due to charged domain walls or strong compensating electric fields [3,4]. Similarly, domain patterns minimize internal mechanical stress through compatible arrangements of spontaneous strain [5]. This results in the well-known 90 • and 180 • domain walls in tetragonal ferroelectrics.
Energy minimization provides a powerful method for understanding how domain patterns form. Since the domain walls are typically less than a few nm in thickness, a sharp interface approximation is attractive [6]. From this approach and from observation of micro-or nano-structure, many domain patterns have been identified [3,5,[7][8][9][10]. The simplest known arrangements are bundles of 90 • and 180 • domains forming laminates of alternating stripes. In thin films or lamellae, the domain bundles are commonly polarized in-plane [11]. However, thicker plates or bulk materials can form energy minimizing patterns with domains polarized in all three axial directions [3,12]. The formation of periodic or nearly periodic patterns of domains is commonly observed; the effect of specimen size and mechanical strain on the formation of such periodic domain patterns has also been widely reported [4,[12][13][14][15]. The spatial distribution of domains affects ferroelectric properties at both fine scale and macroscale [8,11,[16][17][18][19][20][21]. Rödel [21] discussed the dependence of piezoelectric coefficients on the volume fractions of domains, while Weng and Wong [22] showed that the macroscale ferroelectric behaviour is enhanced in certain laminated domain patterns. There is great interest in the domain patterns because of the potential to exploit engineered configurations of domains for energy harvesting, actuation, sensing, and memory devices [23][24][25][26][27].
Significant progress in understanding domains has been made with relatively simple models. However, results based on a sharp interface approach neglect the details of domain walls and their interactions, which become important when the size of domains is comparable to the domain wall width. In such cases, the formation of domain patterns and their stability has not been fully explored. Phase field models, with an order parameter that varies across the domain walls, offer an opportunity to study nanoscale domain patterns with a physically realistic representation of domain wall energy, bulk electromechanical energy, and the effects of external electromechanical loading. The use of phase field models for ferroelectric domains is well-established [28][29][30][31][32][33]. Several studies explore the evolution of polarized domains in thin films [34][35][36][37], under the effect of electromechanical loads. These studies explore the local effects of substrate strain on domain shapes [36], domain wall movement [34] and the coercive field [37]. In ferroelectric-bulk, the 3-dimensional effect of the electromechanical boundary conditions influences the formation of equilibrium domain patterns. 3D periodic laminates with microstructural features such as curved domain walls [5,15] and ribbon-shaped domains [9] have been experimentally observed. These 3D microstructural features evolve under external electromechanical loads and influence the nanoscale behaviour of ferroelectrics [9]. Exploring the spatiotemporal evolution of periodic polarization patterns in a theoretical framework would indicate the complexities in pattern formation and provide initial steps towards nanoscale experimentation. Recently, we used the same phase field model as in this work to explore nanoscale domain patterns [14]. Starting from the patterns predicted by energy minimization with sharp interfaces [5], it was found that several of the patterns dissolve into simpler, more stable arrangements such as single domains or alternating stripes. For computational speed, the calculations were mainly carried out using 2D models that limit the freedom of the system to relax into fully 3D domain patterns. Extending the study to consider periodic volumes large enough to support domain patterns makes heavy demands on computation, and the calculations reported in the present work typically took days to weeks each, running as parallelized code using 32 cores on an Intel Xeon CPU at 3.2 GHz.
The goal of this work is to establish and explore two kinds of periodic boundary conditions for scalar and vectorial fields at the bound of representative volume elements (RVE). This is motivated by the fact that periodic boundary conditions generally restrict the solution space. Thus, any RVE simulation may benefit from this discussion since our model demonstrates that differing periodic boundary conditions naturally lead to an enriched variety of periodic domain patterns. The principle of the reverse periodic boundary condition is an enhanced assembly rule, where we superpose translation and rotation.
We study a range of periodic cell sizes and then focus on a size just sufficient to allow complex domain patterns to form. Such complex patterns could be useful in memory applications where arrays of writable domains and skyrmion-like topologies [38] are of current interest. Periodic patterns may be artificially enforced by imposing a template or may arise naturally through minimization of energy. Nature, of course, has available any periodic cell size or none. The exact reason why nearly periodic patterns of domains commonly form is not known. It may be that similar conditions of stress or cooling rate during the phase transition favour the nucleation of particular domains that then adjust into a nearly periodic form to minimize energy. More likely, a small region of pattern first forms and then becomes the template for surrounding regions which copy the pattern. We show that cell sizes of at least tens of domain wall widths are necessary in order to support domain patterns other than simple stripes or single domains.
Further, we explore the effect of strained states on the formation of domain patterns. The study is motivated by a classification of low-energy periodic domain patterns given by Tsou et al. [5] in which several families of periodic domain patterns were identified. We employ a 3-D phase field model [39,40], which is calibrated for barium titanate and reproduces the transition from an initial paraelectric state to the tetragonal ferroelectric state at the Curie point, whereupon domain patterns evolve to minimize internal energy. We do not force any particular pattern by imposing a template or by initial conditions: instead, the simulations take a randomly perturbed initial state and relax this towards equilibrium. Since it is impractical to simulate 3-dimensional regions much larger than a few tens of nm in size, we make use of periodic boundary conditions, which allow the assembly of simulated volumes into larger periodic cells. Additionally, the simulations are carried out with controls on the average strain, introduced as various side conditions. This allows for the effect of a remanent ferroelastic strain, appearing at the phase transition from the paraelectric to the ferroelectric state.
The paper is organized as follows: the underlying model is given in Sect. 2, followed by a discussion of periodic boundary conditions in Sect. 3. The results of our simulations are presented and discussed in Sect. 4, followed by conclusions in Sect. 5.

The material model
The model follows the Landau-Devonshire theory of ferroelectrics [41], the work of Fried and Gurtin [42,43], and the subsequent work of Landis [39,40]. Fuller details concerning the normalization of physical fields, phase field modelling, and finite element technique are given in [25,44,45]. The evolution of domain patterns follows the time dependent Ginzburg-Landau equation [39], with polarization P as the order parameter: where is the Helmholtz free energy of the system and β is a polarization viscosity introduced to represent the dissipation associated with domain wall motion. Domains of homogeneous electric polarization evolve during the simulations, and the polarization rate on the right hand side of Eq. (1) vanishes in the final state of equilibrium. The Helmholtz free energy includes contributions due to mechanical strain ε, electric polarization P, and electric field E: = mech (ε , P) + well (P) + grad (Grad[P]) + elec (P , E). ( The free energy function used includes six energy wells corresponding to the six polar variants of the tetragonal ferroelectric. The curvature of the wells is calibrated to reproduce the elastic, piezoelectric, and dielectric behaviour of the individual domains, and the gradient term reproduces 90 • and 180 • domain walls with widths and energies consistent with other studies [39,46,47]. The model does not explicitly include flexoelectric effects, though gradients in strain are indirectly coupled to polarization through the gradient terms in the model. Simulations are initiated with P ≈ 0, ε ≡ 0 and E ≡ 0 throughout the simulated region. Strain controls, where used, are switched on after the first iteration of the Newton-Raphson solution scheme. The overall internal energy = B c dV reduces as the simulation proceeds, eventually stabilizing at an equilibrium value. Since the initial energy 0 = B c 0 dV is zero, the internal energy of the model becomes negative during the phase transition. This follows Devonshire's approach [1,2] wherein the free energy of BaTiO 3 is zero in the cubic state. The minimum possible value for is given for the tetragonal crystal structure with homogeneous polarization P 0 in a state free of stress or electric field. Then, the model yields 0 = −0.5345 E 0 P 0 , where P 0 = 0.26 Cm −2 , is the remanent polarization of the relaxed domain and E 0 = 2.182 × 10 7 Vm −1 is the coercive electric field of the single domain (which is much greater than the coercive field of polydomain material). The results given in Sect. 4 are normalized using 0 , E 0 and P 0 .
The model simulates a region of material B c within which electric displacement D is defined by where κ 0 = 8.854 × 10 −12 Vm C −1 is the permittivity of free space, and electric field E is derived from a scalar electric potential φ as The model fulfils Maxwell's electrostatic equations where ρ is the free charge density. Since the material is a ceramic insulator, we assume ρ = 0. This, in combination with Eqs. 3 and 6, penalizes divergent polarization fields, such as 'head-to-head' or 'tail-to-tail' domains. At surfaces with outward normal direction n and surface charge density q, the electric displacement satisfies Equation (5) Since we are not concerned with extrinsic forces, the balance of momentum reduces to while at surfaces with traction t, Periodic boundary conditions on the pseudo-surface, ∂B c , of unit cells will be defined in Sect. 3, regulating the values of P, φ and u on ∂B c . The majority of simulations represent cubic volumes of material, modelled with up to 20 3 = 8000 triquadratic finite elements of approximately 2 nm element size. This yields of the order of 10 6 degrees of freedom. Extra degrees of freedom are introduced at domain walls, as they form, by enriching elements with quartic shape functions wherever the gradients of P exceed a critical magnitude; details of the enrichment method have been given by Muench and Krauß [44,45]. At the initial stage of the simulation, random values of P of magnitude 0.01P 0 are imposed on the system. This gives a perturbed near-zero polarization field corresponding to a state close to the maximum potential energy, representing the cubic state of the crystal. The subsequent evolution of domain patterns is due to the transformation from paraelectric to ferroelectric material at room temperature. Unless otherwise specified, the results show equilibrium states withṖ = 0. Equilibrium was checked at the end of each simulation by increasing the time increment t of the implicit time integration scheme by a factor of 10 6 relative to the time steps used during the simulation, and checking for any change in polarization.

Periodic boundary conditions
Let us consider a nanoscale cubic unit cell B c taken from the interior of a bulk single crystal with periodic domain pattern. Cartesian coordinates x i with orthogonal basis vectors e i aligned to the edges of B c are used. The origin x 0 lies at the centre of B c , as shown in Fig. 1. Dimensions are given by L c such that coordinates c . The surface ∂B c of the unit cell has six faces with surface normal vectors The Each vector r I inherits the parameters of the face ∂B I c as cartesian coordinates and is perpendicular to n I as shown in Fig. 1b.
These periodic boundary conditions enforce zero average strain and electric field within the cell. In Sect. 3.3, we will introduce a side condition to allow for nonzero average strain. By then, Eq. (15) ensures that unit cells are geometrically compatible. Equation (13) also yields solenoidal polarization when averaged over B c : Nevertheless, Eq. (16) does not imply that the polarization is locally divergence-free within B c .

Reverse periodic boundary condition
We next consider an assembly of cubes B c to form an extended unit cell with internal symmetry. This can be achieved by translation along e i combined with a 180 • rotation R I about e i . The periodic conditions then read The 180 • rotation tensors R I are defined by where 11 is the identity tensor and ⊗ yields the dyadic tensor product. With the restriction in Eq. (20), the electric polarization satisfies q ≡ 0 on ∂B c when cubes B c are assembled after duplication, translation along e i via cell length L c , and rotation by R I . This reverse periodic boundary condition can be applied to one, two, or three faces simultaneously. In the present work, the reverse condition is applied on a single pair of opposite faces, with standard periodic boundary conditions employed on the other faces.
As with the standard periodic boundary condition, the reverse periodic boundary condition yields solenoidal polarization averaged over B c . Similarly, this boundary condition implies that the mean strain and electric field vanish over the extended unit cell. Thus, the assembly of cells B c with reverse periodic boundary condition does not require any macroscopic electric field to stabilize it. Figure 2 illustrates the effect of the standard and reverse boundary conditions. The reverse periodic boundary condition encourages the formation of additional 180 • domain walls and hence more complex domain patterns.
In summary, the reverse periodic boundary condition is applied to accelerate the search for complex domain patterns at smaller RVE's. It enforces a double cell symmetry condition, without increasing the computational volume. This boundary condition encourages periodic patterns with oppositely polarized domains to form, in order to maintain an overall electric field neutrality.

Side conditions for strained states
Next, we introduce mechanical side conditions to apply four different representative strained states. We use these representative strained states to explore the general trends in temporal evolution of polarization patterns and to explore the complexities of 3D-domain arrangements.
When BaTiO 3 transforms freely from the cubic phase to the tetragonal phase, the transformation results in a remanent strain ε r , which is non-isochoric. Three variants of ε r are observed, with the tetragonal c axis aligned to e 1 , e 2 and e 3 , respectively. The cubic lattice constant a cub is replaced by tetragonal lattice constants c tet and a tet , respectively, giving remanent strain components Thus, the three strain variants are In each case, the c axis aligns with the polarization P in the minimum energy state. For BaTiO 3 , set a cub = 4.000 Å, c tet = 4, 0328 Å, and a tet = 3.9892 Å, based on experimental data and theoretical results presented in [48]. Then, ε c = 0.0082 and ε a = −0.0027 giving 2 ε a + ε c = 0.0028, which indicates that the unit cell undergoes a non-isochoric process if unconstrained. Thus, it is desirable for the simulations to allow an overall straining of B c , corresponding to displacement of the boundaries ∂B c . The strain-free boundary conditions of Eqs. (15) and (19) would suppress an average remanent strain, resulting in a stressed state. To enable a state closer to stress-free conditions, we introduce scalar parameters λ i , i=1..3, tracking approximate volume fractions of each of the three strain variants, as indicated by the local polarization direction: Using the volume fractions λ i and the remanent strain variants from Eq. (22), define an overall remanent strain Note thatε r is not the actual strain or strain average of the simulation, but instead is a notional remanent strain that is consistent with the current polarization state of the simulation. Since it results from an average over B c , ε r is not a function of position. The overall volume change due to average remanent strainε r is Now consider imposing some state of average strain of the formε ∈ Diag. This could be achieved by replacing the displacement boundary conditions, Eqs. (15) and (19), with for the standard periodic boundary condition, and for the reverse periodic boundary condition. Equivalently, we may keep Eqs. (15) and (19), and replace (8) by the side condition taking care to adjust the definition of Helmholtz free energy in Eq. (2) accordingly. This latter approach is used in the present work. In order to control the average strain during transformation from cubic to tetragonal in differing ways, the three diagonal elements ofε can be defined using parameters f 1 , f 2 , and f 3 by analogy with Eq. (24), such thatε The parametrization with f i aids the formulation of different types of constraints or strain controls. For example, setting f 1 = 1 and f 2 = f 3 = 0 recovers the strain state of a stress-free single domain with c axis aligned to e 1 . By contrast, setting all f i = 0, i = 1 . . . 3, produces a state free of average strain, while setting all f i = 1/3 produces an equiaxed volume dilation that matches the dilation of the cubic-tetragonal phase transformation. Simulations were conducted with four specific side conditions, as follows: First, the average strain-free case, where is denoted Isochoric Ferroelectric Phase Transition (IFPT). Thus, the volume of B c is held constant during the simulation. A negative spherical stress is likely to arise to restrainε r if the material evolves polarized domains such that volume fractions λ i > 0. By contrast, a low stress phase transition may be promoted by Here,ε adapts continuously to matchε r ; this condition will be denoted Adaptive Spherical Stress Free (ASSF). Note that ASSF may generally yield unequal strains along the three axes, but tracks any volume dilation accompanying the transformation, so that typically tr[ε] = (2 ε a + ε c ) = 0.0028 at equilibrium. A third condition enforces equal average straining along each axis while matching the expected volume dilation due to transformation: This will be denoted as the Isotropic Spherical Stress Free (ISSF) case. This side condition tends to eliminate any average spherical stress. Finally, a common arrangement for materials bound by a substrate, or subjected to a uniaxial mechanical loading, is that the strain state retains two equal principal strains. Combining this with the restriction to match the volumetric dilation due to the phase transformation gives: The particular case when b = 1/2 we denote Transversely Isotropic Deformation (TID). Note that all four cases allow for locally inhomogeneous strain within B c . Further, the volume fractions λ i need not generally coincide with imposed strain parameters f i , except for the case of ASSF, where there is continuous tracking to match these parameters. The strained states might bear similarities to experimentally imposed mechanical strains. However, the spirit of the current work is to indicate the general trends for complexity and types of patterns that form in bulk ferroelectrics.

Results and discussion
Domain topologies resulting from the use of periodic boundary conditions in 3-dimensions, combined with each of the four strain conditions (IFPT, ASSF, ISSF, TID), are presented in this section. It is instructive first to consider the evolution of domains within a single simulation, then to explore the effect of the periodic cell size in the range 9-41nm on the topologies that may form. Finally, a cell size large enough to generate complex domain patterns (41nm) is selected and the domain topologies resulting from each set of boundary conditions are described.

Evolution of a 2-dimensional domain topology
In this example, a cube B c of size L c = 33 nm is simulated using 16 × 16 × 16 elements. Standard periodic boundary conditions are used on the opposite face pairs (∂B 1 c , ∂B 4 c ) and (∂B 2 c , ∂B 5 c ). Meanwhile, the reverse periodic boundary condition is applied to the face pair (∂B 3 c ,∂B 6 c ). The adaptive spherical stress-free condition (ASSF) on average strain is applied. Figure 3 shows a series of states of the simulation as it evolves towards equilibrium; each image shows four copies of the simulated volume B c to illustrate the periodic continuation in the e 1 and e 3 directions. The reverse periodic boundary condition encourages (but does not enforce) the formation of a 180 • domain wall normal to the e 3 direction.
In the early stages of the simulation, the local values of |P| are much less than P 0 , but the energetic driving force for the formation of polarized regions rapidly produces domains. These initially nucleate in a random, high energy pattern, and gradually resolve sharply defined 90 • and 180 • domain walls. As can be seen from Fig. 3c-e, the domain pattern is complex as it evolves towards equilibrium, but eventually reaches a simpler, lower energy state.
The final equilibrium state shows a herringbone-type pattern that was also identified by Tsou et al. [5], using a sharp interface model, as a member of the family {1314} of compatible rank-2 domain patterns. In this case, the final pattern is essentially 2-dimensional, with the polarization vector always lying in the e 1 − e 3 plane; this was not forced, but rather was a natural outcome of the relaxation towards a minimum energy state. The same pattern was noted in earlier work using a phase field model, but in a 2-dimensional simulation with a 40nm periodic cell size [14]. However, in that study a fixed average strain was imposed to stabilize the pattern. Here, the adaptive strain condition allows the topology to form and adopt an average strain state that minimizes energy.

Cell size variation with standard periodic boundary conditions
In this section, the effect of periodic cell size by varying the cell size L c from 9 to 41 nm is under investigation. Below about 9 nm size, the cell is too small to support multiple domains (domain wall widths in the model are in the region of 2-3 nm. Beyond 41 nm the computation becomes excessively cumbersome due to the scaling of the number of degrees of freedom in the problem. From our previous work on scale effects in polarization patterns [13], we expect that increasing the cell size will enable polarization patterns with multiple domain walls to form. However, in the current work we stop at 41nm cell size because it is sufficient for complex domain patterns to form within practical computation time. For the smaller cell sizes (9-25 nm), elements of size 1 nm were employed. Cell sizes of 33 and 41 nm were simulated with 16 3 and 20 3 elements, respectively. The standard periodic boundary condition is used as a representative boundary condition to explore the effect of cell size variation. We expect a similar trend in the results with reverse periodic boundary conditions. Several runs of the model were performed with different random initial polarization perturbations; these runs can reach different equilibrium states. The results in Fig. 4 show the final domain topologies with lowest overall Helmholtz free energy. Additionally, the side condition defining the average strain was varied to explore each of the four cases defined in Sect. 3.3.
As shown in Fig. 4, simple 90 • domain stripes arise in most of the simulations. However, the TID constraint, which imposes a strain condition close to the free strain of the single domain, strongly encourages a monodomain to form. Generally, smaller cell sizes result in simpler domain patterns. Only the 41 nm cell is large enough that more complex domain patterns can form if this reduces the energy. Thus, for the subsequent study of domain topologies, the 41 nm cell size is employed.
Next, let us define volume fractions α i , β i , i = 1, 2, 3, to assist in identifying the presence of 90 • and 180 • domains: where the integrals are defined over the entire periodic cell volume V c (note that this may contain multiple, transformed copies of cube B c ). In Eq. (34), V p gives the volume fraction of the periodic cell that is polarized into domains. At equilibrium this is normally close to 100%. The values of α i and β i give, respectively, the fractional net polarization along each axis and the fraction of domains aligned with each axis, normalized in each case by the total volume fraction of domains. For convenience, we sort these volume fractions such that  Table 1 shows the normalized Helmholtz free energy along with α and β as percentage values for the equilibrium states, computed with 41 nm cubes.
The normalized overall free Helmholtz energy = B c /L 3 c dV is always negative, as the zero datum is the high energy, cubic, state. However, the ratio / 0 , where 0 = −0.5345 E 0 P 0 corresponds to the energy  Table 1. Imposing a strained state (IFPT, ISSF, TID) often does not lead to the single domain solution as can be seen from the α, β values; even when a monodomain does result, the energy is higher than that of the stress-free state. Allowing the periodic cell to adapt its strain as domains develop (ASSF) results in the lowest energy solutions, which may be a simple monodomain, but can develop a domain pattern if the random starting conditions nucleate one. The isochoric conditions (IFPT) result in the highest vales of free energy, as the periodic cell is rigidly constrained to match that of the cubic state.

Domain patterns simulated with standard periodic boundary conditions
From now on, we only consider results from simulations of 41 nm cells. Figure 5 shows the domain pattern for the isochoric ferroelectric phase transition (IFPT) and using the standard periodic boundary conditions. Since these conditions strongly constrain the formation of domains aligned with any one preferred direction, it is expected that a fully 3-dimensional pattern of domains with some regions polarized parallel to each of the three directions n I (I = 1 . . . 3) is likely to form. Tsou et al. [5] classified the simplest compatible laminate of domains that has all three tetragonal strain variants present as a rank-2 laminate in a family {1325} that has a mixture of 90 • and 180 • domains. This pattern is also the outcome of the phase field simulation: Fig. 5b-d shows cross sections of the pattern that identify the polarized directions of the domains.
Close examination reveals that the laminate does not have perfectly planar interfaces between domains, but instead the intersections of domains form complicated angular shapes that result in vortex-like regions at the triple junctions where domains polarized along the ± e 3 and − e 2 directions meet. There has been great interest in such vortex centres [11,38,49] because of their potential use in applications such as memory storage, and their interesting topological properties as skyrmions. The simulations here suggest that it may be possible to generate arrays of vortices by suitably constraining crystals during the cubic-tetragonal phase transition. Changing the mechanical side condition to the isotropic spherical stress-free case (ISSF) keeps the average strain spherical (equal in all directions) but reduces the energy by allowing volumetric expansion to match the phase transformation, see Table 1. The result shown in Fig. 6 is once again a rank-2 domain pattern that was identified in the study of Tsou et al. [5]. In contrast to the IFPT solution in Fig. 5, it does not contain 180 • domain walls. However, this pattern, in the family named {1315}, contains an electrically incompatible boundary where domains polarized in the − e 1 direction meet domains polarized in the − e 2 direction. As can be seen in Fig. 6, this frustrated boundary adopts a curved shape to reduce energy. Thus, the equilibrium state is once again an imperfect laminate that would not be expected from consideration of the constrained, sharp interface theory alone. We note that domain arrangements observed on the cross sections of the multi-cell assembly in Figs. 5 and 6, bear resemblances to polarization patterns in thin films, simulated using phase field methods [51]. For example, the vortex-like regions at triple junctions and the curved domain walls are predicted to form on the free surface of a ferroelectric thin-film geometry. The transversely isotropic deformation condition (TID) enforces a strain close to that of the unstressed single domain, and the simulation then produces a monodomain equilibrium state as shown in Fig. 4. However, the adapted spherical stress-free strain condition (ASSF), which continuously adapts the strain to be consistent with the average polarization, allows other patterns to form. The monodomain is a possible solution, but in our simulations the ASSF condition more commonly resulted in periodic 90 • domains, as shown in Fig. 7. The TID and ISSF cases have lower energy than that of IFPT, but the adaptive case (ASSF) gives the lowest energy of all, since it relaxes the average strained state as far as possible while retaining periodicity.
In summary, the model tends to form more complex domain patterns, including some of the vortex structures and ribbon-like domains shown in Figs. 5 and 6, if the mechanical side condition of periodic boundary conditions constrains the ferroelectric phase transition. The reason for this is that the constraint forces the material to produce a mixture of domains to match the imposed strain. In experimental work, the constraint can be imposed by the use of epitaxial growth with lattice mismatch. This method has been used in strontium titanate, for example, to provide biaxial tensile stresses that induce incipient ferroelectricity [52]. By contrast, our work suggests that lattice mismatch or stress gradients in the sample could be used in thin films of ordinary ferroelectrics to promote complex domain structures, or to encourage specific patterning of domains [53]. For example, a stress gradient generated by an atomic force microscopy tip was used to mechanically write polarization patterns on a barium titanate thin film [54]. We show exaggerated deformations in Fig. 8a to illustrate that, although the volumetric strain is held at zero, there is nevertheless local straining of the periodic cell. The pattern has roughly equal volumes of domains with their c axis aligned parallel to the e 1 , e 2 and e 3 directions, as indicated by β = [39,33,28]. Vortex centres are evident at the junctions of the e 3 and ± e 2 domains in Fig. 9d. Further simulations using the isotropic spherical stress-free side condition (ISSF) yield similarly complex patterns, see Fig. 10. These domain patterns have many 180 • domain walls, along with ribbon-like domains [50] that thread through the structure. In both cases, five of the six polarization variants are present, which can be understood as follows: the strain condition enforces the presence of some domains with c axis aligned to each of the three coordinate axes. However, polarization must balance in both the ± e 1 and ± e 2 directions due to the rotation of the cube during assembly, so all four variants polarized in the (x 1 , x 2 ) plane must be present. Meanwhile, since some domains with c axis aligned to the e 3 are expected at least five variants are likely to be present overall. From the methods employed by Tsou et al. [5,6], it is immediately evident that the pattern formed cannot be any laminate of rank less than 3. The patterns shown appear similar to rank-3 laminates, but with imperfect alignment of domains and curved domain walls. The equilibrium domain patterns in Figs. 9 and 10 both exhibit centres of vorticity and ribbon-like domains, with five of the six polarization variants present. Close inspection reveals that they are related variants. This is reflected in their similar α and β values (see Table 1).
The result of applying the adapted spherical stress-free case (ASSF) and reverse boundary condition on the 41 nm cube is the same as that shown in Fig. 3g for the 33 nm cube. Thus, there is no need for an additional figure.
Finally, the result of applying the transversely isotropic deformation case (TID), along with the reverse periodic boundary condition, is shown in Fig. 11. As expected, the formation of a 180 • domain wall is enforced. However, the result is as simple as possible: we observe alternating 180 • domain stripes, a well-known domain pattern in barium titanate. While the TID constraint with reverse periodic boundary condition also admits a single domain equilibrium state, it appears that the alternating stripes is a more likely end state because of the high probability that, during the phase transition, different regions in the periodic cell may nucleate domains aligned in different directions. As the simulation proceeds towards equilibrium, these regions grow and stabilize. There is then an energy barrier to be overcome if one region is to be eliminated to allow the pattern to become a monodomain. As with the standard periodic boundary condition, the equilibrium states computed with the reverse periodic boundary condition show a progression in free energy values, reducing in the order IFPT→TID→ISSF→ASSF, consistent with the gradual reduction of constraint from the most restrictive (ε = 0) to the least restrictive, adaptive condition.

Impact of an outer electric field
Finally, we demonstrate the effect of an outer electric field on the stability of a 3D ribbon-like periodic domain pattern solution. As representative example, let us choose the polarization pattern modelled with isochoric constraint (IFPT) and standard periodic boundary condition. The electric field is applied along the + e 3 direction. In Fig. 12, the volume fractions of the six polarization variants in the RVE are plotted as coloured lines. The corresponding arrangements of the domains are shown as subset images at significant stages.
In the initial stages, up to E = 0.08E 0 the smaller domains with polarization in the − e 3 direction shrink in size and disappear. However, the domains with polarization in e 1 and − e 2 direction remain such that the ribbon-like domain pattern is stable up to E = 0.64E 0 , beyond which a striped domain pattern with e 1 and e 3 domains forms. On increasing the electric field up to E = E 0 , the domains with polarization along e 1 direction remain in the polarization pattern. We interpret this polarization pattern as arising from the competing isochoric mechanical side condition (IFPT) and the applied electric field. At E = E 0 , the IFPT constraint stabilizes the e 1 -domain in the presence of opposing external field. However, at E = 1.21E 0 , the stripe pattern dissolves to form a monodomain with net polarization along e 3 -axis. The net polarization in the domain pattern decreases during the transition from the stripe pattern to a monodomain. We associate this decrease in net polarization to result from the strong IFPT mechanical constraint. The mechanical strains modelled in the transverse direction to the applied electric field, force the monodomain away from its relaxed domain state.

Conclusions
In this work, we have applied multiple periodic boundary conditions for the simulation of ferroelectric domain patterns on the nanoscale. For the well-known standard conditions, our solutions tend to evolve patterns with 90 • domain walls only, if the RVE is allowed to adopt the inherent volumetric growth resulting from the materials phase transition from the cubic to the tetragonal crystal state. This result has been observed for three important variants (ASSF, ISSF, TID) of many more possible mechanical side conditions. However, we have also found a mechanical side condition enforcing 180 • domain walls. It restricts the RVE concerning its isochoric deformation, which we have denoted the isochoric ferroelectric phase transition (IFPT).
Since many experimental works do not restrict the isochoric deformation of the crystal during its ferroelectric phase transition, and since they have found ferroelectric patterns with 180 • domain walls, we have proposed the reverse periodic boundary condition, which naturally enforces 180 • domain walls within the RVE by an enhanced continuation rule of the cells. This extends the capabilities of our simulations and reproduces some domain patterns that are well known from prior theoretical studies and from observations of BaTiO 3 .
The reverse periodic boundary condition might be of interest for similar simulations in multiphysics; we just discuss one specific application here. But even for pure mechanical problems it widens the solution space of the displacement field. We refer to Fig. 8a, where the exaggerated deformation of the equilibrium state is deliberately shown. Obviously, displacement fields exist such that the assembly of cells requires the superposition of translation and rotation. To the best knowledge of the authors, the reverse periodic boundary condition is novel in the context of RVE and ferroelectric simulation models.
For both kinds of periodic boundary conditions, the simulations also exhibited interesting features such as frustrated domain walls, curved ribbon-like domains, and centres of vorticity. Some of these features can be of use in ferroelectric devices, and the simulations suggest that, by controlling strain during the cubic-ferroelectric transition, it may be possible to generate arrays of domains with these exotic features.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.