Towards Inﬁnite Tilings with Symmetric Boundaries

: Large-time coarsening and the associated scaling and statistically self-similar properties are used to construct inﬁnite tilings. This is realized using a Cahn–Hilliard equation and special boundaries on each tile. Within a compromise between computational effort and the goal to reduce recurrences, an inﬁnite tiling has been created and software which zooms in and out evolve forward and backward in time as well as traverse the inﬁnite tiling horizontally and vertically. We also analyze the scaling behavior and the statistically self-similar properties and describe the numerical approach, which is based on ﬁnite elements and an energy-stable time discretization.


Introduction
At first-order phase transitions, coexisting macroscopic domains of different phases emerge from small fluctuations of a homogeneous phase. Late stages of this process are often dominated by the motion of the interfaces separating the domains. Considering this large-time coarsening behavior, i.e., the growth of a characteristic length scale l(t) as t → ∞ determines important characteristics of the dynamics and led to the identification of several universality classes of domain growth. We are here concerned with conserved order parameters, for which the expected behavior is l(t) ∼ t 1/3 , which results from the scale invariance of the Mullins-Sekerka system x → λx, t → λ 3 t. Rigorous results exist for an upper bound for l(t), stating that microstructures cannot coarsen faster than the similarity rate [1]. As there are non-generic configurations, e.g., stripe domains with zero curvature which are stable, lower bounds cannot be expected within a deterministic framework. Besides this scaling law, solutions with random initial data are also believed to be statistically self-similar in this large-time regime. Numerical studies based on a Cahn-Hilliard equation and related coarse-grained theories indicate that the approach of the large-time regime with the statistically self-similar structures might be very slow [2]. To explore these regimes numerically thus requires large length and time scales, which limits the accessible sample size. We are here interested in these statistically self-similar structures, which have been used for various art and design projects, e.g., [3]. Here, we would like to explore very large, in principle infinite, samples. To tackle such a system we consider, instead of one huge simulation, many moderately sized domains with different initial data, and require the boundaries to match. If appropriately done, this will allow construction of large (infinite) tilings which are statistically self-similar. With a random arrangement of finitely many computed structures, the impression of an infinite tiling with no recurrence could be achieved. For this impression, the boundary conditions at the computational domains are crucial. They are described in detail in Section 4 together with the finite-element approach to solve the Cahn-Hilliard equation. In Section 2 we show various results, among other things a computer program which allows navigation through space and time of an infinitely extended structure. We further discuss improvements and outline possible applications. In Section 3, we discuss scaling and self-similar properties.

Results
As the underlying model for phase transitions with a conserved order parameter we consider a Cahn-Hilliard equation see Section 4 for details. First, we consider coarsening in a rectangular cuboid using standard boundary conditions n · ∇φ = n · ∇µ = 0 on ∂Ω. Figure 1 shows snapshots of the results within the large-time regime, visualized in various ways. The interface area is minimized, and the structure thus coarsens in time. To analyze this in a statistical manner requires either a larger domain or more samples. Our idea is to combine both by using samples which are distinct from each other but fit together to form a large and extendable structure. The boundary conditions in the current setting enforce the level lines of the interface to be perpendicular to ∂Ω. They in addition do not fit to each other and thus do not allow combination of different samples. To overcome these limitations, we consider a smaller domain, again a rectangular cuboid and the boundary conditions introduced in Section 4 which specify the values of φ and the normal flux ∇φ · n such that opposite sides match. The first approach only has two distinct boundaries, N = 2. Figure    Even if the inner structure is unique for each tile, the boundaries in xand y-direction are the same and recurrences are visible. To improve on this issue, we consider a second approach, which is less flexible in terms of arrangement of samples but minimizes possible recurrences. As a compromise of computational cost and visual impression we consider tilings with N = 10 different boundary conditions. This improves the impression of a tiling with no recurrence since systematic recurrence in a row, a column, or diagonally can be avoided by careful assembly of the tiles. Repetitions do not only appear less frequently but also in a pattern that is much less obvious. To see this recurrence without knowing the construction process and explicitly searching for them is almost impossible. We consider two different internal realizations each, leading to M = 10. Within the proposed pattern determined by the boundary conditions the inner realization are randomly chosen to construct an infinite tiling, where recurrences are almost invisible.
A software is developed to visualize the infinite structure. We consider visualizations with five projected level lines of the interface. The software allows zooming in and out, evolve forward and backward in time as well as traverse the infinite tiling horizontally and vertically. Figure 3 shows some screenshots, starting from an early time instant and a low zoom factor (a), going to a late-time instant of this setting (b), zooming into the structure (c), evolving along a trajectory in space and time, which keeps the interface area constant (d), and going back to the initial state (a). Videos of the journey through space and time are provided in SI as Videos S8 and S9. As the proposed approach is in principle not restricted to rectangular cuboidal domains various possibilities for applications can be imagined. Besides wallpaper design, they range from fashion design with individualized clothes to camouflage patterns of automotive prototypes. Here, we highlight a more entertaining application, a Rubik's cube which always fits, but has 24 different fields; see Figure 4.

Discussion
We now explore the scaling properties and the statistically self-similar behavior of the constructed infinite tiling. The theoretically scaling behavior l(t) ∼ t 1/3 of a characteristic length scale is tested by computing the interface area of each tile over time. The interface is thereby represented as the level-surface φ = 0. In addition, we also compute the length of the interface of the level line at z = 0.25, which is used for visualization. Figure 5 shows the results over time, which are averaged over all individual realizations of tiles. The results lead to a scaling exponent which is below the upper bound of 1/3. The slope is not constant, but on average equal for the two measures and approximately 0.26. There are different reasons for this lower value, either the structure has not reached the late-time behavior for which the theoretical scaling behavior is expected, or the considered domain with the zero-flux boundary conditions on top and bottom favor parallel structures and thus prevent coarsening. The limitation of our approach, to combine several tiles and to compute them separately, should also be mentioned. This approach only allows the consideration of coarsening up to a length scale of the size of a tile. For larger times, the approach is no longer valid. However, even if the theoretical scaling law could not be shown computationally, statistical self-similarity still might be possible. Statistical self-similarity can already be expected from a randomly chosen sample. We consider the middle slice at an early time instant, its coarsening and subsequent zooming-out in Figure 6. Instead of the interface line the two phases are rendered in black and white. When the coarse structure is zoomed out to a degree where the interface-length matches that of the earlier timestep, both structures are visually similar (see Figure 6a,c). To quantify this, we compute for each row and column of pixels the distance between two interfaces. This is done for all samples and plotted for different times in Figure 7.  Even if the theoretically predicted scaling law l(t) ∼ t 1/3 could not be computationally shown, the large-scale simulations, which run for each tile on a high-performance computer, reproduce the predicted statistical self-similarity. The huge structure, which results as an arrangement of individual tiles, would not have been possible to simulate on the available hardware. The approach fulfills two goals, it provides enough statistics to analyze scaling and statistical self-similarity and it allows the generation of aesthetically appealing tilings with almost invisible recurrences, which can be infinitely extended.

Materials and Methods
The Cahn-Hilliard equation [4] is a fourth order partial differential equation resulting as a H −1 gradient flow of a Ginzburg-Landau energy where Ω is a bounded domain, a positive parameter determining the width of the diffuse interface, γ the surface energy, here considered as a positive constant, and B(φ) = 1 4 (1 − φ 2 ) 2 a double well potential. The resulting equation reads and converges for → 0 to the Mullins-Sekerka problem [5]; see Ref. [6]. Various numerical approaches have been proposed to solve the equation efficiently. We consider a convexity splitting approach, e.g., [7][8][9][10][11]. The idea is to split the double well potential B(φ) = B c (φ) − B e (φ), such that both parts are convex and to consider the time discretization as with discrete time derivative d τ φ n+1 = (φ n+1 − φ n )/τ n . The resulting scheme is unconditionally energy stable, unconditionally solvable and converges optimally in the energy norm [9]. To solve the above systems, we consider a linearization of We further consider adaptive mesh refinement according to criteria related to the position of the diffuse interface, here ∇φ, to ensure a resolution of approximately five grid points across the interface and a coarser mesh elsewhere; see Figure 8. The resulting linear system is solved in parallel using a block-preconditioner, see [12,13], and the iterative solver FGMRES. All problems are implemented in the adaptive finite-element toolbox AMDiS [14,15]. The considered parameters are = 0.01 and γ = 1.0 and as computational domain rectangular cuboid Ω = (0, L)x(0, L)x(0, l) with l = 0.2 and L = 2.0 for the larger and L = 1.0 for the smaller domain and boundaries Γ top , Γ bottom , Γ 0 and Γ 1 = 4 i=1 Γ 1,i , see Figure 9. The number of grid points on the larger domain reduces from approx. 2.95 million at the beginning to approx. 1.95 million at the final configuration. As initial condition we consider white noise around the mean value φ = 0. At Γ top and Γ bottom we specify zero-flux boundary conditions for φ and µ. We first consider the larger domain with zero-flux boundary conditions for φ and µ also on Γ 0 . In this setting the finite-element formulation in each time step reads: with V h = {v ∈ C 0 (Ω)|v |T ∈ P 1 (T)∀T ∈ T } the space of piecewise linear Lagrange elements and triangulation T . In each time step we extract φ n+1 1 = φ n+1 and n · ∇φ n+1 1 = n · ∇φ n+1 along the inner boundary Γ 1 . These data are going to be used in the subsequent computations in the smaller domain as boundary conditions on Γ 1 . The finite-element formulation now reads in each time step: For different initial data this leads to different solutions with common boundary conditions. However, to construct tilings, they also must match, which is not yet guaranteed. To fulfill this requirement, we proceed in two different ways. The first approach considers only one computation on the larger domain and uses the extracted values and fluxes φ n+1 1 and n · ∇φ n+1 1 at Γ 1 only from two sides Γ 1,1 and Γ 1,2 (N = 2) and specifies them also on the opposite sides for the computations on the smaller domain. For M different initial conditions this generates M individual samples which match with each other at the boundaries if translated by the domain size in xor y-direction. This leads to a very flexible arrangement of the samples but has the drawback of frequent recurrence at the boundaries.
The second approach also begins with a computation on the larger domain with boundary Γ 0 but subsequent computations are performed on intermediate domains that extend to the bounds of Γ 0 in directions where a fresh structure is desired at the boundary and are restricted to the bounds of the smaller domain Γ 1 in directions where the structure is to be continued from an already existing neighbor tile by using its stored values and fluxes φ 1 and n · ∇φ 1 at Γ 1 . This requires small modifications of the finite-element formulation in Equations (7) and (8) using only parts of Γ 1 instead of the whole inner boundary. When enough samples are computed to define all N boundary sides the remaining samples are computed on Γ 1 with all sides fixed by boundary conditions from earlier computations.
The most simple setup meeting our design-goal of non-obvious recurrence of boundaries requires five different tiles A 0 , B 0 , C 0 , D 0 and E 0 which can be assembled into a row that matches the same row displaced by a few tiles above and below. Then, five rows of five tiles each form a square which can be continued in all directions indefinitely (see Figure 10c). This setup allows for ten unique boundary sides s 1 through s 10 (see Figure 11). Inside our 5 × 5 square of tiles recurrences do not occur in a row, a column, or diagonally, which would not be possible with a smaller number of unique tiles. With only two or three different tiles, repetitions would have to occur at least diagonally (see Figure 10a) and with four different tiles it is only possible to build a unique 4 × 2 rectangle without diagonal repetitions (see Figure 10b).
On the macroscopic scale our pattern still has obvious repetitions since we need to continue the same square of 5 × 5 tiles in all directions to form the infinite tiling. To remedy this problem, we generate variants of our five initial tiles A 1 , B 1 , C 1 , D 1 and E 1 which exactly copy the boundaries of their respective prototype with index zero but differ on the inside. Now, in our infinite tiling each tile of a specific prototype is replaced randomly with either the original 0-variant or the new 1-variant of that type. With only two variants per type we already allow for more than 33 million (2 25 ) distinct squares of 5 × 5 tiles. For our interactive visualization software, we settled with these M = 10 tiles due to memory-constraints. However, for printed realizations such as wallpapers, further variants could be computed, making it even harder to spot recurrences.
The initial five tiles must be computed consecutively since their boundaries depend on each other. Starting with A 0 we have no constraints yet and thus the simulation yields four fresh boundary sides. The simulation domain Γ 0 exceeds the region of interest Γ 1 in all four connecting directions because we need to avoid the level lines always being perpendicular to the boundary of Γ 1 (Figure 11a).  As a next step we simulate a temporary tile B tmp to prepare the run for B 0 . Obviously, tile B 0 has to match A 0 's right-hand side boundary (s 2 in Figure 11b) at its own left-hand side. However, this is not the only constraint. Our 5 × 5 layout requires that B 0 meets A 0 also in its top-right corner (see Figure 10c). Thus, in our temporary step besides fixing the left boundary we also fix the top boundary with a mirrored version s m 3 of the data from A 0 's bottom side. The mirroring is required because the top-right corner of B 0 must conform to the bottom-left corner of A 0 . The simulation domain for B tmp exceeds the region of interest only in two directions: bottom and left, where we obtain data for two fresh boundary sides (Figure 11b). Now we are ready to compute our second tile B 0 . We fix the right, bottom and left boundaries with data from our previous two simulations and this time leave the top-side unconstrained to obtain a fresh boundary that is only fixed at the two corners where it will match tile A 0 (Figure 11c).
For tile C 0 we require another preliminary run C tmp to fix the bottom-left corner where it meets A 0 . We use the mirrored top-side s m 1 of A 0 at C tmp 's bottom to account for that and fix the left and top sides where C 0 shares sides with B 0 and A 0 . Only at the left-hand side we obtain a fresh set of boundary-data (Figure 11d). For C 0 we release the bottom-constraint of C tmp and obtain another fresh boundary (Figure 11e). The last tile with a fresh boundary is D 0 . All sides, except the right-hand side one, are already constrained by previous computations (Figure 11f). For E 0 all four sides are fully predetermined (Figure 11g). We need to simulate it only for its interior.
All variant-tiles A 1 through E 1 (and further ones if desired) can now be computed in parallel since their boundaries are already known. They are all restricted to Γ 1 and do not require extraction of boundary-values and -fluxes, hence are computationally less expensive than the initial runs.
Funding: This research was funded by DFG grant number SFB/TRR96 A7. We acknowledge computing resources provided by the Jülich Supercomputing Center under grant number HDR06.