Tuning the Stability of a Model Quasicrystal and its Approximants with a Periodic Substrate

Quasicrystals and their periodic approximants are complex phases, which have by now been observed in many metallic alloys, soft matter systems, and particle simulations. In recent experiments of thin-film perovskites on solid substrates, the type of complex phase was found to change depending on thermodynamic conditions and the type of substrate used. Here, we investigate the effect of a substrate on the relative stability of a two-dimensional model quasicrystal and its approximants. Our numerical methods are molecular dynamics simulations and free energy calculations that take into account phason flips explicitly. For weak substrates, we observe an incommensurate-commensurate transition, in which a continuous series of QC approximants locks into a discrete number of approximants. For stronger substrates, an enhancement of the stability of the dodecagonal quasicrystal and a variants of square lattices were found. All phenomena can be explained by the interplay of the model system with the substrate. Our results demonstrate that designing novel complex periodic and quasiperiodic structures by choice of suitable substrates is a promising strategy.


Introduction
The discovery of structures with point group symmetries forbidden in periodic crystals 1 created a new field of crystallography, the field of quasicrystals 2,3 (QCs).Bulk, three-dimensional QCs became popular as candidates for materials with functional properties such as thermal insulators 4 , coating materials 5,6 , photonic crystals 7 and superconductors 8 .The first atomic QC in two dimensions was found only recently in the thin-film perovskite BaTiO 3 on a Pt substrate 9 .This discovery inspired follow-up works of thin-film perovskite QCs as well as their approximants at different chemical compositions [10][11][12][13][14][15] .Despite significant efforts, how these perovskite QCs form and why they are stable remain mostly open questions.Answering these questions is of interest to simulation studies but difficult to address with quantum mechanical methods at the atomic scale because the lack of periodicity and slow phason relaxation in QCs requires large systems over long simulation times 16,17 .Coarse-grained model systems are advantageous because they capture the essential physics and offer a route towards elucidating relative thermodynamic stability of the involved phases and their phase transformation pathways.Two-dimensional particle systems in the presence of a substrate  Institute for Multiscale Simulation, Friedrich-Alexander-Universität Erlangen-Nürnberg, 91058 Erlangen, Germany.† E-mail: michael.engel@fau.deexhibit intricate structure formation phenomena caused by the competition of two length scales.These length scales are the size of the unit cell of the unperturbed particle system and the periodicity of the substrate.The competition of the length scales can introduce novel commensurate and incommensurate phases 18,19 .Incommensurability also plays a role in nanotribology, for example Aubry transitions in systems of colloids sliding on a substrate [20][21][22][23] .Three factors are relevant: substrate potential depth, substrate periodicity, and the crystallographic symmetries of the particle system and the substrate 24 .
The investigation of thin QC films on a substrate can be divided into two classes of systems.The first class comprises systems, which are templated to have QC order by use of a QC substrate.Such QC films have been realized at the nanoscale as atomic adlayers [25][26][27][28] .And at the colloidal scale, QC laser fields can be finely tuned to induce a gradual transition from an unperturbed system to a frozen QC [29][30][31][32][33] .The second class comprises systems with intrinsic QC order.That is, such systems are QCs even in the absence of the substrate.It is suspected that thin film perovskite QCs fall into the latter class 9 .Both classes of systems can show interesting phase behavior as the substrate is tuned.
There exist many simulation models for 2D QCs.Prominent examples are binary mixtures of particles interacting with the Lennard-Jones potential 34,35 , one-component systems of particles with interaction potential minima at two length scales [36][37][38] , and anisotropic particles with patchy interactions 39,40 .Here, we study the Lennard-Jones-Gauss (LJG) potential, which stabilizes a number of one-component QCs and approximants already in the absence of a substrate 37,41 .We use this coarse-grained model to investigate the possibility of targeting specific approximants by tuning substrate parameters.As we will see, the LJG system can behave like a thin film QC in both the first and the second class, depending on the choice of thermodynamic conditions and the substrate.

Coarse-grained model
As in our previous work 41 , we investigate a coarse-grained twodimensional model system, which is known to form different types of (quasi-)periodic phases 37 .In this model, identical particles interact isotropically with the LJG potential given by . ( The potential parameters are fixed in this work at ϵ LJG = 1.8, r LJG = 1.42, and σ 2 LJG = 0.042 to stabilize a dodecagonal QC at intermediate temperatures 41 .The phase diagram of the LJG potential at these parameters contains a number of phases described by tilings made from squares, triangles, and pentagons.The sequence of phases can be described according to the underlying symmetry change 41 p3m1 0.24 where the symbols are the wallpaper groups of the phases and the numbers above the arrows indicates the critical temperature of the phase transitions.It has been shown that a relatively simple periodic approximant (wallpaper group p3m1) is the energetic ground state stable at T = 0.This simple approximant transitions via a series of increasingly more complex hexagonal approximants (wallpaper group p6m) continuously to a dodecagonal QC (wallpaper group p12m) in the temperature range 0.24 < T < 0.35.Towards slightly higher temperature, at T = 0.36, the QC converts into a square tiling (wallpaper group p4m) before melting at T = 0.41 into the liquid.The liquid has no broken rotational or translational symmetry and is thus represented by the group of all isometries, the Euclidean group E(2).
The phase behavior of the LJG system is only known in the absence of a substrate.Here, we add a substrate by including an external potential.We mimic the effect of the Pt (111) substrate for thin-film perovskite QCs by choosing an external potential in form of a hexagonal lattice by superposition of three plane waves, plane (see assembly diagram below, Fig. 5).A schematic depiction of the model, comprising of the particles (red) interacting with themselves and the periodic substrate is shown in Fig. 1.
To justify the choice of model system for this work, we need to discuss phason modes.Phason modes are dynamic modes characteristic for QCs.They are realized as collective rearrangements of finite clusters of tiles and can be decomposed into a series of elementary particle displacements called phason flips.In many common QCs, like for example the Penrose tiling, phason mode excitation requires only low activation energy and can occur locally with only a small number phason flips.We call phason modes such as those found in the Penrose tiling continuously excitable.In other QCs, for example the square-triangle tiling, phason mode excitation requires creating a pair of defects and propagating each defect along a path rearranging the tiling in the process until the defects annihilate.In the square-triangle tiling this process is called a zipper. 42Zippers are not efficient in relaxing the square triangle-tiling.We call phason modes such as those found in the square-triangle tiling not continuously excitable.Continuously excitable phason modes are found in thin-film perovskite QCs where they involve rhomb tiles with small interior angle 30 • .Because of the presence of continuously excitable phason modes, perovskite QCs can efficiently relax their tiling and transform to and from approximants 10,11 .We chose to study the LJG model system because it stabilizes includes continuously excitable phason modes in the p6m and p12m phases 41 .This means our model system is similar to perovskite QCs in the character of its phason flips and the phase transformations between approximants and the QC.

Tiling construction in hyperspace
The tilings of the dodecagonal QC and its approximants in the LJG model system can be described by the cut-and-project method 43,44 in a four-dimensional hyperspace.We construct tiling models in two steps.
In a first construction step, approximants of the quasicrystalline shield tiling 45 are constructed.The projection matrices on the two-dimensional physical (or parallel) space is with c  = cos(( + 1/ 2)π/ 6) and s  = sin(( + 1/ 2)π/ 6).The projection matrix on the two-dimensional internal (perpendicular) space is with c ′  = cos(π/ 6 + ϕ) and s ′  = sin(π/ 6 + ϕ).The occupation domain for the shield tiling is a dodecagon with twelve vertices constructs the dodecagonal shield tiling for ϕ = 0 and q/ p periodic approximants for ϕ ̸ = 0.In this notation, approximants are labeled by two integers, q and p.Their ratio ranges from q/ p = 2 for the simple periodic approximant with wallpaper group p3m1 to q/ p = ( 3 + 1)/ 2 for the quasicrystal with wallpaper group p12m.This means the simple periodic approximant is identified as the 2/ 1 approximant.All approximants with wallpaper group p6m can be constructed in the range ( 3 + 1)/ 2 < q/ p < 2. The approximants have rectangular unit cells with lattice constants b  = b y 3 and b y = (t + 3)p/ 2. The approximants of the shield tiling constructed so far are built from three basic tiles: square, triangle, and shield 45 .
In a second construction step, the basic tiles square, triangle, and shield are decorated by placing one particle in the center of the square tiles, one particle in the center of the triangle tiles, and three particles as a regular triangle in the center of the shield tiles.The result are the tilings of interest, built from three larger tile types: square, triangle, and slightly deformed pentagon.Examples of these resultant tilings are found in Fig. 1 and below in Fig. 6.

Simulations and free energy calculations
We perform simulations and free energy calculations using a combination of molecular dynamics (MD) and Monte Carlo (MC) to accelerate relaxation and exploration of phase space.MD is used to sample phonon modes and configurational entropy.MC is used to sample phason modes 46 .From the van Hove autocorrelation function 38 , we observe that all phason flips readily occur over a single flip distance near r = 1 and that most particles can perform phason flips.This demonstrates that phason modes are continuously excitable via local phason flips and do not require significant activation energy.We use the HOOMD-blue molecular dynamics simulation package 47,48 as a basis for our simulations and further extend it to include MC moves.For MD, the integration timestep is 0.01.The NVT ensemble with Nosé-Hoover thermostat is used and open boundary conditions are applied.The Boltzmann constant is set to 1.For MC, we attempt a phason flip in MC by choosing a particle at random and a flip vector uniformly distributed in the ring 0.9 < r < 1.1.The MC move is accepted according to the Metropolis acceptance criterion.The measured acceptance rate of MC phason flips is about 10%.Periodic boundary conditions are employed to minimize finite size effects.We now discuss the two different simulation modes used in this work.
In the first simulation mode, used for generating a diagram summarizing assembly behavior as a function of model parameters, called an assembly diagram, and for studying phase transitions, we perform MD simulations without MC moves.We simulate 10000 particles in a quadratic simulation box of box edge length 300.Particles are initially placed in a central circle with number density 0.3.The particles crystallize over time into a roughly circular, compact cluster, which corresponds to a solidgas coexistence.We chose this setup to effectively obtain open boundaries for the crystalline cluster.Such open boundaries are essential for efficient phason relaxation because phason relaxations are generally suppressed by the topological constraints of periodic boundary conditions.
In the second simulation mode, used for calculation of Helmholtz free energies, we combine MD and MC.5000 MD sweeps, i.e., moving each particle 5000 times, are alternated with 5000 MC steps, i.e., attempting to perform a phason flip for 5000 randomly chosen particles.System size varies for different approximants and the QC between 4000 and 12000 particles.In contrast to the first simulation mode, the simulation box for approximants is now completely filled to prevent a solid-gas coexistence and to avoid phase transitions between approximants.For the QC, we simulate with open boundaries and only consider particles sufficiently fare from the surface.We first relax constructed approximants in an NPT simulation and then perform free energy integration in the NVT ensemble.In this way, the number density of simulations in the second simulation mode is guaranteed to agree with that of simulations in the first simulation mode.
We use the Frenkel-Ladd (FL) method 49 to compute free energies.FL performs a thermodynamic integration over a parameter λ.It interpolates between the system of interest at λ = 0 and an Einstein crystal at λ = 1.In an Einstein crystal, each particle is connected to a reference position by a harmonic spring.The integration is performed as where ΔF is measured from simulation.Because the reference free energy of the Einstein crystal, F Ein , is known, we can calculate in this way absolute free energies.
It is important to check that the system remains in equilibrium throughout the FL integration.We perform such a check by calculating free energy in forward direction (increasing λ) and in backward direction (decreasing λ) and by checking for hysteresis.Absence of hysteresis is an indication that free energy calculation is accurate.A test FL integration for the 2/ 1 approximant is shown in Fig. 2. Forward and backward curves for ΔF agree, except for a small difference near λ = 0.At small λ, particles diffuse strongly through the system, which slows down equilibration.We use our test FL integration to adjust the numerical parameter to ensure that hysteresis is sufficient small.For numerical integration in production calculations, we use Legendre-Gauss quadrature with 20 λ points.Each λ point is simulated for 25000 MD sweeps.The spring constant of the Einstein crystal is Λ E = 100.This value guarantees similar vibrational dynamics in the Einstein crystal and the system of interest as indicated by an approximately flat integrand 50 .
We validate the FL implementation by comparison to reported free energies 41 for the continuous transition p3m1 → p12m in the absence of a substrate at ϵ = 0. We calculate free energies F n (T) for six constructed approximants q/ p = n/ 10 with n ∈ [15, . . ., 20] and a well-equilibrated QC.We follow the continuous transition from the 2/ 1 approximant to the quasicrystal via the superstructure wave vector 41 This wave vector is calculated from the distance between firstorder satellite peaks and main diffraction peaks in reciprocal space 41 .The appearance of satellite peaks is a consequence of the presence of structural modulations in the approximants.We define the stability temperature of the approximant n/ 10 as the intersection of F n−1 and F n+1 .As an exception, instead of the 14/ 10 approximant, we use a well-equilibrated quasicrystal.
Our implementation of the FL plus phason mode relaxation algorithm shown in Fig. 3a accurately reproduces prior work 41 .In particular, we reproduce the reported linear increase of the wave vector q s as a function of temperature (Fig. 3b).The agreement with prior work validates our free energy calculation with the FL plus phason model relaxation method.

Molecular dynamics simulations
We analyze the crystallization behavior of our system by performing 77 MD simulations at seven temperatures 0.2 ≤ T ≤ 0.5 in steps of 0.05 and for eleven substrate potential depths 0.0 ≤ ϵ ≤ 0.1 in steps of 0.01.Substrate periodicity is fixed at λ = 0.5.This value is inspired by a direct comparison of relevant length scales in our simulations and the perovskite QCs on a Pt substrate 11 .Each simulation is run for 10 8 integration steps and the final simulation snapshot at the end of the run is analyzed.In this subsection, we first perform a general analysis of the orientational order and then conduct a more in-depth analysis of the observed structures by visual inspection and with the help of diffraction patterns.

Local orientational order analysis
We quantify the presence of n-fold local orientational order via the Mermin order parameter 51 defined as where N is the number of particles, N j the number of nearest neighbors of particle j using the cut-off r = 2.5, and θ jk the orientation of the nearest neighbor bond connecting particle j and particle k.The Mermin order parameters for n = 6 and n = 4 are also called the hexatic and quadratic order parameter, respectively.We are specifically interested in these order parameters because we expect to observe a competition between the quadratic order at high temperature T (as previously found in the LJG potential 41 ) and the hexagonal order of the substrate.
The results of the local orientational order analysis are shown in Fig. 4. Indeed, there exists a competition between quadratic and hexagonal local orientational order.Hexagonal local order dominates at high T.Here the liquid phase is expected.In contrast, quadratic local order dominates at intermediate T and high ϵ.Apparently, the substrate does not promote local hexagonal order but instead local quadratic order.This observation will be discussed further below.The gradual variation of local orientational order across the T-ϵ parameter plane indicates that the particles can react smoothly to changes of temperature and the substrate at a local level.We also observe that the presence of thermal fluctuations dominates over crystallographic symmetry in our two-dimensional system, smearing out the boundaries separating ordered phases.This precludes using the Mermin order parameters for automatic calculation of a assembly diagram without relaxation or temporal averaging.

Classification of observed phases
To obtain an assembly diagram, we analyze each final simulation snapshot in direct space by inspection of the particle configuration as well as in reciprocal space via calculation of its diffraction pattern.The latter is important because approximants can exhibit significant phason fluctuations (see discussion below).We   manually compare snapshots to tiling structures with the following wallpaper groups: p3m1 (2/ 1 approximant), p6m (other hexagonal q/ p approximants), p12m (dodecagonal QC), p4m (square lattice), and E(2) (liquid).Manual classification is typically unique and unambiguous.In the presence of multiple tiling structures or in case defects remain in simulation snapshots, a snapshots is classified as belonging to the dominant tiling structure.
We summarize the structure classification in the assembly diagram of Fig. 5.In the absence of a substrate, ϵ = 0, the known phase behavior of Eq. ( 2) is reproduced.Notice that due to low resolution along the temperature axis, the continuous series of q/ p approximants reported previously 41 is not fully resolved.Phase boundaries shift with ϵ according to four trends: First, the dodecagonal QC and the 3/ 2 approximant become stable over broader temperature ranges and shift towards lower temperature.This finding demonstrates that a hexagonal substrate can, in fact, stabilize a dodecagonal quasicrystal despite it having lower symmetry than the quasicrystal.Second, the stability region of the square lattice shifts towards high T with increase of ϵ, and a new structure, called the modulated square lattice (msq), appears.Third, melting temperature increases with ϵ.Fourth, at low T and high ϵ we find a novel chiral quasicrystal (cqc).In cqc, the inversion symmetry of the dodecgonal quasicrystal is broken.Because this manuscript focuses on weak substrate potentials, cqc will be subject of a follow-up work.Three simulation snapshots at parameters marked by asterisks in Fig. 5, the 2/ 1 approximant at (T, ϵ) = (0.25, 0.00), the 3/ 2 approximant at (0.25, 0.02), and the quasicrystal at (0.25, 0.04) are shown in Fig. 6.A fourth simulation snapshot at parameters marked by an asterisk at (0.3, 0.09) corresponding to msq is shown in Fig. 8.
We make two further observations: First, in the absence of a substrate, the phase transformation from the 2/ 1 approximant to the QC was continuous via a series of q/ p approximants.This continuous behavior is not sustained towards finite ϵ.Instead, higher-order approximants are generally suppressed.Specifically, we identify the low-order approximants 2/ 1 and 3/ 2 repeatedly in our snapshots for ϵ > 0 but no other approximants.Second, there remain significant phason mode fluctuations in the 3/ 2 approximant at elevated temperatures.Such phason mode fluctuations are apparent when comparing Fig. 6b-left, which shows a simulation snapshot, to Fig. 6b-right, which shows the ideal tiling of the 3/ 2 approximant obtained after slow relaxation to T = 0 to remove phason fluctuations.Notice the similarity of the diffraction patterns calculated for both snapshots (two halves in Fig. 6b-middle).Apparently, the 3/ 2 approximant behaves like a random tiling in that it permits phason mode fluctuations in thermodynamic equilibrium.This means, we must rely on characterization of average crystallographic order by identification of diffraction spots in the diffraction pattern to distinguish the approximants in the generation of the assembly diagram.
The presence of phason mode fluctuations in an approximant is unexpected and to the best of our knowledge has not been reported.It is unclear whether such behavior is related to the dimensionality of our system and can also occur in three dimensions.We observe random tiling behavior only in sufficiently com- plex approximants.For example, in the 2/ 1 approximant there are no square tiles (green squares in Fig. 6) except if connected to structural defects.This means phason flips cannot occur inside a 2/ 1 approximant but must be initiated from the surface or from a grain boundary.In contrast, the 3/ 2 approximant has many square tiles and thus readily supports phason mode fluctuations.

Incommensurate-commensurate transition
While the unperturbed (i.e., in the absence of a substrate) 2/ 1 approximant transforms into the dodecagonal QC continuously via a series of hexagonal approximants, all approximants except 2/ 1 and 3/ 2 quickly disappear when the substrate is turned on (Fig. 5).Why is this the case?To understand the disappearance of higher-order approximants, we resort to free energy calculations.We calculate free energies for a well-relaxed QC and a number of low-order approximants.Recall that hexagonal approximants are labeled by two integers as q/ p with ( 3 + 1)/ 2 ≈ 1.37 < q/ p < 2. Without loss of generality, we can require q and p to be coprime.Finally, q and p should be small such that the approximants have small unit cells.The simplest such approximants are 2/ 1 = 2, 9/ 5 = 1.8, 5/ 3 ≈ 1.67, 3/ 2 = 1.5, and 7/ 5 = 1.4.We arrange the approximants in increasing order of q/ p ratio from the 2/ 1 approximant, which is the energetic ground state, to the QC, which the equilibrium state at intermediate temperature.QC in correct (that is, in q/ p decreasing) order with temperature.We also confirm the disappearance of all approximants except 2/ 1 and 3/ 2 with increasing ϵ.Note that we did not include in the free energy calculations the square lattice and its modulation and the liquid.Besides this omission, there is good agreement of the phase diagram in Fig. 7 with the assembly diagram obtained from MD in Fig. 5.
Situations where the competition of two length scales dominates phase behavior have long been studied in the field of modulated crystals 52 .Because the two length scales are of different physical or chemical origin, their ratio is not expected to be any special number and typically will be irrational.The requirement of the system to satisfy both length scales simultaneously then causes the emergence of structural complexity in form of modulated crystals.One distinguishes two cases: Incommensurately modulated crystals are aperiodic crystals where both length scales are realized without compromise.The dodecagonal QC in our system can be interpreted as an incommensurately modulated crystal.A related example are Moiré patterns, where not a competition of length scales but the competition of crystallographic orientations and translation generates interference effects.Incommensurately modulated crystals are accompanied by phason modes in thermodynamic equilibrium.In contrast, commensurately modulated crystals are periodic crystals with large unit cells.The lattice constants of commensurately modulated crystals simultaneously approximates integer multiples of both length scales.Commensurately modulated crystals result from a lock-in transition where the irrational ratio of the two length scales locks into a rational approximation.q/ p approximants are examples of commensurately modulated crystals.
Phase diagrams that include modulated crystals are usually studied as a function of length scale ratio and coupling strength.
The first parameter, length scale ratio, is often controlled indirectly by changing the composition or local structure.In our system, length scale ratio is linked to temperature T because temperature affects tile energy and thus the tile type that is dominant.The second parameter, coupling strength, is measured relative to thermal energy.In our system, coupling strength is given by substrate potential depth ϵ.With these identifications, we can compare our phase diagram to past phase diagrams of modulated crystals.The comparison shows that Fig. 7 closely resembles the mean field phase diagram for the 3D Ising model with competing interactions (Figure 24 in Ref. 53 ) and the FVdM model as suggested by Aubry (Figure 25 in Ref. 53 ).Such characteristic phase diagrams of incommensurate-commensurate transitions have been called the Devil's flower phenomenon by Bak 54 .Our results confirm that the increase of substrate strength locks modulated crystals into certain approximants and drives new complex phase behavior as expected from the Devil's flower phenomenon.
The appearance of the Devil's flower in our phase diagram helps us interpret the observed stabilization mechanism and particle dynamics.In the absence of a substrate in the temperature range 0.24 < T < 0.36, the system is incommensurately modulated.Phason modes are unlocked and critical for thermodynamic stabilization because they contribute to free energy.As substrate potential depth increases, phason modes may still be present (see Fig. 6) but gradually disappear.Eventually the system locks into the 2/ 1 and 3/ 2 approximants.To the best of our knowledge, Fig. 7 is the first report of an incommensuratecommensurate transition in any two-dimensional system of freely moving particles.

Square lattice and its variants
The presence of the substrate drives the incommensuratecommensurate transition at low substrate potential depth ϵ < 0.03.For ϵ > 0.03, the square lattice (sq) and its variants appear dominantly in the phase diagram (Fig. 5).For intermediate substrate potentials, a novel modulated square lattice (msq) is observed (Fig. 8a,b).For high substrate potential, ϵ ≫ 1.0, above the range covered in Fig. 5, a coexistence of grains with rectangular lattice (rec) and hexagonal lattice is eventually observed (Fig. 8c,d).All of these phenomena are yet other examples of commensurate-incommensurate transitions induced by the competition of the length scales imparted by the LJG interaction potential with the unit cell of the substrate.
Msq and rec are typically polycrystalline with individual grains locking into one of three distinct orientations.This locking is a consequence of the six-fold symmetry of the substrate.Diffraction patterns confirm the six-fold symmetry.In the diffraction pattern of msq (inset of Fig. 8a), isolated peaks (circled) are present for each grain.Besides the locking of grain orientations, individual grains are only weakly affected by the substrate (see modulations discussed below).In the diffraction pattern of rec (inset of Fig. 8c), some peak triplets merge into single peaks (indicated by black hexagons in Fig. 8a,c).The merging is caused by a stretching of squares into rectangles to improve compatibility with the substrate.This means four-fold symmetry is already lost in individual grain of rec and the system is closer to six-fold symmetry than msq.While small patches of triangles are sometimes already present in rec (central part of Fig. 8c labelled '4'), only very high ϵ transforms the system completely into a hexagonal crystal (not shown).
A zoom-in on grain '2' of rec is shown in Fig. 8d.We plot the substrate potential in the background to visualize particle locations (red) relative to maxima (yellow) and minima (green) of the potential.Potential maxima form a hexagonal lattice separated by channels of low potential that trace the nearest-neighbor bond network of a honeycomb lattice.The stretching of squares into rectangles can now be understood.Recall the nearest-neighbor distance is r NN ≈ 1.0, and the substrate periodicity is chosen as λ = 0.5.A lattice constant   = 2λ allows particles to retain their preferred nearest-neighbor distance without penalty.And a slight adjustment of the lattice constant to  y = 4/ 3  ≈ 1.15  ensures commensurability with the substrate also in the y-direction.
A zoom-in on grain '2' of msq is shown in Fig. 8b.Just like before, a lattice constant   = 2λ allows particles to retain their preferred nearest-neighbor distance without penalty.The y-direction is again more complicated.A lattice constant  y =   now guarantees four-fold symmetry in the diffraction pattern but can do so only in average and by allowing positional modulations along the y-axis.A close inspection reveals the presence of these modulation in form of weak satellite peaks in the diffraction pattern (inset of Fig. 8b).The separation of satellite peaks to the central peaks corresponds to a superstructure periodicity of  s ≈ 0.20.Analysis of other snapshots suggests that  s is not unique but a function of T and ϵ.We also observe that particles are unusually mobile in channels along the y-direction (purple in Fig. 8b).
This high mobility can be directly associated with a phason mode.All evidence combined demonstrates that msq is an incommensurately modulated crystal.

Conclusions
Research in this paper was originally motivated by the discovery of approximants in the thin-film perovskites (Ba, Sc)TiO 3 on solid metal substrates.Importantly, experiments demonstrated that the type of approximant changes as a function of substrate interaction and periodicity.Our findings reproduce and explain aspects of these experimental observations.Our findings also highlight a number of new phenomena, not yet seen in experiment or simulation.
The new phenomena we observed, and that we believe could be general consequences of the competition between a QC and a periodic substrate, are: 1. Existence of incommensurate-commensurate phase transitions.In these phase transitions, the QC locks into one or several approximants.We observed the stabilization of the 3/ 2 approximant via a lock-in transition.
2. Approximants with phason fluctuations.We observed that the 3/ 2 approximant resembles a random tiling rather than an ideal crystal despite having a well-defined periodicity as indicated by the presence of a reciprocal lattice.
3. Enhanced stability of the QC.Counter-intuitively, the substrate enhanced the thermodynamic stability of the dodecagonal quasicrystal for some parameter values.
4. Incommensurately modulated crystals, where the modulation is a consequence of the competition between two length scales.We observed the modulated square lattice.
5. Novel structures commensurate with the substrate that are not approximants.We observed the rectangular lattice.
While observation 1 is similar to that in thin-film perovskites and observation 3 could explain the existence of perovskite QCs, the other observations have not yet been reported in the experimental system.It will be interesting to search for them.Generally, we see a trend that the competition between the inherent symmetry of the model system, i.e., the crystal structure the system prefers to form in the absence of a substrate, and the imposed symmetry of the substrate can occur in multiple steps involving both commensurate and incommensurate phases.We find in our system, as we increase the strength of the substrate potential, at T = 0.25 (see Fig. 5) a complex sequence of phase transformation alternating multiple times between commensurate and incommensurate: from (i) 2/1 approximant (crystal) to (ii) 3/2 approximant (commensurate) to (iii) quasicrystal (incommensurate) to (iv) modulated square lattice (incommensurate) to (v) rectangle lattice (commensurate) to (vi) hexagonal lattice (substrate).
Future work should further expand this research by testing more parameters or conducting similar studies in more realistic systems.An important parameter that should be varied, but was

Fig. 1
Fig. 1 Schematic depiction of the model system of this work.Particles (red) interacting with themselves via the LJG potential and with an external potential simulating the substrate.The particles form a dodecagonal quasicrystal.Bonds connect nearest neighbors and form three types of tiles: triangles, squares, and pentagons.Tiles are deformed due to vibrational motion of the particles at an elevated temperature.a, Top view.b, Side view at an angled perspective.

32 Fig. 2
Fig.2Frenkel-Ladd integration for the 2/ 1 approximant at T = 0.25.We measure the integrand ΔF in simulation by time-averaging at each λ point in both forward direction and backward direction.The inset shows a zoom-in to small λ values.

Fig. 3 aFig. 4
Fig. 3 a, Helmholtz free energies F n (T) for a well-equilibrated quasicrystal and six approximants n/ 10, n = 15, . . ., 20.Standard errors are estimated from eight independent FL calculations.Purple plus markers indicate the stability temperatures of approximants.Data shows the forward branch of the free energy calculation.b, Validation of the free energy calculation with the FL plus phason model relaxation method.The superstructure vector grows linearly with temperature from the 2/ 1 approximant to the QC.Error bars on the temperature axis connect the FL calculation for the forward and the backward calculation.Stability temperatures for approximants agree with reference values 41 within the error bars.

Fig. 7 Fig. 6 7 TFig. 7 a
Fig. 6 Snapshots from MD simulations at T = 0.25 and at three substrate potential depths ϵ.Nearest neighbor bonds are shown only, with the particles omitted.Triangle tiles are shown in red, square tiles in green.a, 2/ 1 approximant at ϵ = 0.00.b, 3/ 2 approximant at ϵ = 0.02.c, QC at ϵ = 0.04.Snapshots are shown in direct space obtained from simulation (left column) and from representative ideal tilings (right column).Diffraction patterns (middle column) are calculated and compared for both snapshots (left and right halves).The reciprocal lattice of diffraction spots in reciprocal space is mapped onto characteristic lattice vectors in direct space as indicated by blue lines in a, b.Notice the presence of significant phason mode fluctuations in the simulation snapshot of the 3/ 2 approximant (b, left).Despite these phason mode fluctuations, its diffraction image (b, left half of middle) is highly similar to the diffraction image of the ideal tiling of the 3/ 2 approximant without any phason mode fluctuations (b, right half of middle).

Fig. 8
Fig. 8 Variations of the square lattice at higher substrate potential depth ϵ and intermediate temperature T. a,b, Modulated square lattice (msq) at ϵ = 0.09 and T = 0.3 for 10000 particles.c,d, Coexistence of small grains with rectangular lattice (rec) and hexagonal lattice at ϵ = 1.0 and T = 0.4 for 1024 particles.a,c, Particle configurations showing nearest neighbor bonds only, with the particles omitted.Triangle tiles are shown in red, square tiles in green.There coexist grains in three orientations (labelled '1', '2', '3') and a region dominated by triangles (labelled '4' in (d)).Diffraction pattern are shown as bottom right insets with diffraction peaks encircled and colored according to the grain orientation they correspond to.b,d, Zoom ins on the grain '2'.Particles are shown as red disks overlaid on the substrate potential.In (b), light blue color indicates the channels the particles move in.A diffraction pattern is shown as top left inset.