Phase Field Crystal model for particles with n-fold rotational symmetry in two dimensions

We introduce a Phase Field Crystal (PFC) model for particles with n-fold rotational symmetry in two dimensions. Our approach is based on a free energy functional that depends on the reduced one-particle density, the strength of the orientation, and the direction of the orientation, where all these order parameters depend on the position. The functional is constructed such that for particles with axial symmetry (i. e. n = 2) the PFC model for liquid crystals as introduced by H. L\"owen [J. Phys.: Condens. Matter 22, 364105 (2010)] is recovered. We discuss the stability of the functional and explore phases that occur for $1 \leq n \leq 6$. In addition to isotropic, nematic, stripe, and triangular order, we also observe cluster crystals with square, rhombic, honeycomb, and even quasicrystalline symmetry. The n-fold symmetry of the particles corresponds to the one that can be realized for colloids with symmetrically arranged patches. We explain how both, repulsive as well as attractive patches, are described in our model.

Sketch of the preferred orientation of neighboring particles with (a,b) n = 4-fold rotational symmetry showing cases that are typical for all even n and (c,d) n = 5-fold rotational symmetry, which is typical for odd n. In (a,c) neighboring particles prefer to possess the same orientation, while in (b,d) they prefer alternating orientations, which we realize by a modulated alignment interaction. As a consequence, the cases in (a,d) correspond to patchy particles with attractive patches, while in (b,c) the patches are repulsive.

Introduction
The formation of patterns in particulate systems is a long-standing topic [1]. A mean field approach to describe the formation of complex equilibrium phases is the so-called Phase Field Crystal (PFC) model [2,3] where a free energy expansion in a density-like field and its gradient is considered in a way similar as in the well-known approaches by Swift and Hohenberg [4], by Alexander and McTague [5], or by Lifshitz and Petrich [6]. Non-trivial phases are stabilized in PFC models in various ways, e. g., by using more than one length scale [7,8] (similar as in the Lifshitz-Petrich model [6,9,10]) or by introducing a competition with an incommensurate external potential [11,12]. In another approach anisotropic particles are considered. These would break the flow of the text, as they do not introduce an orientational field. Instead they directly couple the director field to the density field. By introducing an orientational field and couplings between the orientational and the density-like field, PFC models for particles with polar [13] or axial [14][15][16] symmetry have been modeled. Here we want to generalize these models for particles with other n-fold rotational symmetries. Note, there are other PFC approaches to describe systems composed of particles with a certain rotational symmetry [17,18]. However, in these approaches the orientation field and the density field cannot vary independently and as a consequence some phases like plastic crystals and oriented crystals cannot be distinguished.
In this work we want to consider both, attractive and repulsive patches. In case of attractive patches the patches of neighboring particles tend to point towards each other while repulsive patches tend to be oriented away from neighboring patches. As a result neighboring particles either have the same orientation or they are rotated by an angle π/n as illustrated in figure 1.
Note that the particles in this work do not possess any hard core. Therefore, large overlaps can occur and the ordering that we report corresponds to those of so-called cluster crystals [33][34][35][36]. Cluster crystals occur naturally in PFC models or related approaches [10] and can form periodic as well as aperiodic structures [10]. In the conclusions in section 4 we will discuss how our approach might be modified in order to describe particles with hard cores, i. e. particles that cannot overlap significantly.
The article is organized as follows: In section 2 we introduce the PFC model for n-fold particles and explain how we determine the stable phases numerically. In section 3 we first present an overview of the phases that occur for various rotational symmetries before we discuss these phases in more detail. Finally, we conclude in section 4.

Model
We generalize the PFC models that have been introduced for particles with 1-fold [13] or 2-fold [14][15][16] symmetry to particles with n-fold symmetry. We will consider both, attractive and repulsive patches. Concerning the notation, we will follow the model presented in [15], which is shortly outlined in the next subsection.

Short summary of the PFC model for 2-fold symmetry
In [15] Achim et al. have proposed and studied a PFC model for apolar liquid crystals, i. e. particles of 2-fold rotational symmetry. The used free energy functional model had previously been derived by Löwen in [14] from classical density functional theory. The free energy F[ψ, U ] is given as a functional of the density-like field ψ and the orientation field U . The density-like field ψ gives the deviation from the mean density and the orientational field is defined in terms of the nematic order parameter ψ 2 and nematic director field ϕ as (1) Hence the modulus of U represents the intensity of orientational ordering, while the complex phase encodes the direction. In slightly modified notation, the free energy is with the complex derivative operator The values of the parameters B l through F can be linked to the temperature, mean density, and the interaction of the modeled particles. In the scope of PFC modeling they can be considered free parameters. The equilibrium phases are found by minimizing the free energy with respect to the fields ψ and U . In the minimization, the mean of the density-like field is conserved, whereas the orientation field is treated as a non-conserved field. Therefore the pseudodynamical equations read with the Lagrange multiplier λ(t), keeping the mean of ψ constant.

New generalization for n-fold symmetry
We now want to generalize the free energy functional of (2) to n-fold symmetry. First, the definition of the orientation field from (1) is generalized to In this way the interpretation of U is compatible with n-fold symmetry, since the complex phase of U is 2π/n-periodic in the director field ϕ. This means U performs a full turn in the complex plane when the particle orientation rotates by 2π/n. Note that the term proportional to F (F -term) in (2) is the only contribution to the free energy that depends on the complex phase of U . Since all terms in F[ψ, U ] need to be invariant under local rotations of the coordinate system, the F -term needs to be adjusted to the new definition of U . The operator κ transforms under a local rotation of the coordinate system by an angle α as κ = κe −iα . Therefore, U κ n is the simplest term that is linear in U and that respects local n-fold rotational symmetry. Hence the F -term of (2) is generalized to A modification of the E-term proves to be necessary for numerical stability as we will explain in more detail in our discussion of the stability presented in section 3.1. Hence we replace the operator −∇ 2 in the E-term by f ∇ 2 . Moreover, this generalization permits modulated alignment that will be used to reflect both repulsive as well as attractive patches as will be explained in section 2.4.
With the E-term and F -term modified, the free energy becomes This is the new free energy functional for particles with n-fold rotational symmetry that we study in this article. When keeping f ∇ 2 = −∇ 2 and specializing to n = 2, our model reduces to the one from [15]. Moreover, the generalization proposed here is in accordance with a model for polar liquid crystals, i. e. particles of 1-fold rotational symmetry (n = 1), presented by Wittkowski et al. in [13]: For example, the F -term for n = 1 can be identified with the term proportional to B 1 in [13]. Table 1 provides a detailed comparison of all terms. Note that although [13] considers polar particles of 1-fold symmetry, contributions of 2-fold (nematic) symmetry are taken into account in [13], by introducing fields for both, polar and nematic order parameter. This is reasonable, since the 2-fold symmetry is compatible with the underlying 1-fold symmetry. The 2fold contributions can thus be seen as an extension upon the lowest-order orientational ordering. Yet, the number of parameters is drastically increased due to the multitude of cross terms. We narrow down the number of terms in the free energy functional by restriction to a single orientation field, representing the fundamental contribution of n-fold symmetry. Table 1. Comparison of terms occurring in the excess free energy for polar (1fold) particles [13], including nematic (2-fold) contributions, with the respective terms for purely 2-fold [14] and 1-fold symmetry. The single terms are referenced by their coefficients in the notation of the respective publication. Terms that cannot occur in a model are marked with a dash (-) and deviations are specified. The terms of [13], proportional to E 1 through G 7 , are unparalleled in [14] and this work. Moreover the terms in [15] are identical to those in [14], up to a rescaling, which merges B and C of [14] into Bx of [15]. n = 2 [14] n = 1&2 [13] n = 1 [this work]

Minimization scheme and numerical details
Motivated by the numerical investigation of systems with axial symmetry by Achim et al. [15] we implement a similar combination of explicit Euler integration in direct space for the contributions from F dir and implicit Euler integration in reciprocal space for the contributions from F rec . Conveniently, the terms of F[ψ, U ], that are quadratic or bilinear in ψ or U , can be expressed in reciprocal space, where spatial derivation reduces to multiplication with wave vector components.
To be specific, the explicit Euler integration in direct space is given by and the implicit Euler integration in reciprocal space leads to an evolution according to  where F denotes the Fourier transformation and L is a matrix that contains coefficients stemming from the functional derivatives of F rec . The terms in the matrix L will be shortly discussed in the stability analysis in section 3.1. After the implicit integration, the Lagrange multiplier is applied, setting the appropriate value for the k = 0component of Fψ . The backwards transformation to direct space completes the timestep.
We apply periodic boundary conditions and discretize ψ and U on a grid of 512 × 512 or 1024 × 1024 points. The intrinsic lengthscale is set by the B x -term to k = 1, i. e. a peak-peak distance of 2π in direct space. The system size spans 10 to 14 peak-peak distances, such that single peaks, as well as structures much larger than the typical unit cell are resolved. We perform several computations for each point in phase space, starting from random noise in both fields. When the resulting free energies differ, which happens when defects or metastable structures appear, we start more computations with appropriately patterned initial fields. The structure that yields the lowest free energy is assumed to be the equilibrium phase.

Modulated or unmodulated alignment to implement attractive or repulsive patches
Here we explain how the function f ∇ 2 (or f −k 2 in reciprocal space) is chosen such that either the alignment or the anti-alignment of neighboring particles is preferred. Note that the differences in alignment of neighboring particles is used to model repulsive or attractive patches as sketched and explained in Fig. 1. The asymptotics of f −k 2 , i. e. the behavior at large wave vectors, are dictated by the stability requirement which will be discussed in section 3.1.
In case neighboring particles prefer to align, we use a monotonic function namely with m = max {2n − 4, 2}. Note that for n = 2 the term used in [15] is obtained.
In case the anti-alignment of neighboring particles is preferred, the function f −k 2 should be non-monotonic leading to a negative free energy contribution for an orientation that alternates on the length given by the nearest neighbor distance. As anti-aligned particles experience a rotation of π/n from one particle to its neighbor, the period of U amounts to 2 particle distances. This corresponds to a wavelength of k = 1 2 . However, we aim not to introduce a second lengthscale via the preferred wavelength of alignment. Thus we choose an extended in f −k 2 , covering k = 1 2 and k = 1. Furthermore, the asymptotics should be similar as for the alignment case. Based on these requirements we employ a continuous piecewise-defined term namely We call this case the modulated alignment case. Since f (0) = 0, the parameter D has the same meaning for modulated and unmodulated alignment and is not shifted by E.

Stability requirement
To study the stability we have a closer look at the implicit integration step (12) that corresponds to the linear contribution of the integration and is governed by a matrix with We now want to choose f ∇ 2 such that the integration is stable. Note that in general the choice that is used in [14,15] is not suitable: For example, if n ≥ 4 and the E-term was given by f ∇ 2 = −∇ 2 , the computations would exhibit a numerical instability for many (significant) regions in physical parameter space.
In figure 2 we plot the matrix elements of (15) as functions of k. All of the matrix elements have poles at certain values of k. These poles derive from a root of the determinant of 1−∆tL, which enters into (15) as denominator. The position of the poles depends on the physical parameters of the free energy functional and also on the timestep ∆t of Euler integration. For sufficiently small ∆t, the poles fall outside the range of numerically accessible wave vectors. Hence ∆t can be chosen small enough to avoid numerical instability caused by the poles themselves.
Yet a different cause of instability remains, independent of ∆t: In certain directions in reciprocal space the diagonal elements that correspond to the curves labeled by 5 in figure 2 exceed 1 for all wave vectors up to the pole. These matrix elements let the high-frequency modes grow, while the terms in F dir of higher order in ψ and U either overcompensate or fail to compensate for this growth.
To guarantee stability we want the respective matrix elements to fall below 1, i. e.
This relation can be used to compute the maximum F , for which the calculations will be stable, given the timestep ∆t, the maximum length of contributing wave vectors k, and the parameters of the free energy functional. In addition, the relation (20) permits us to pinpoint a condition on the asymptotics of the E-term, necessary for numerical stability at arbitrary non-zero F . Let f −k 2 be of order O(k m ) for large wave vectors k 1. Then (20) is expressed as In order to ensure that this relation holds for a finite F and arbitrarily large k, m ≥ 2n − 4 is required. Conversely the function in the E-term has to behave asymptotically like or like a term of higher order in k. Only then the right hand side expression of (21) does not tend to zero for large k.
Note that qualitatively the order of the term f −k 2 seems not to matter as long as the integration is stable. For example, for liquid crystals (n = 2) we have found the same phases as in [15] even if we use f −k 2 = k 4 instead of f −k 2 = k 2 .

Overview of the phase behavior
Due to the number of parameters and the rich zoo of complex phases that can occur, systematic studies of the whole parameter space are difficult. Our goal in this work is not to present complete phase diagrams. In this and the following subsections we want to present typical phases that can occur in systems for various rotational symmetries. Concerning the parameters we have focused on the regions in parameter space where non-trivial phases are expected as discussed e. g., in the work by Achim et al. on 2-fold particles [15].
In figure 3 we show an overview of the phases that we have found in systems with 1 ≤ n ≤ 6 both, for alignment between neighboring particles (left column in figure 3) as well as modulated alignments (right column). While the figures shown in figure 3 not necessarily are complete phase diagrams as we cannot be sure whether there are additional intermediate phases (see also discussion in section 3.4), the figures demonstrate the typical phases that occur in the small and large limits of D for almost all rotational symmetries namely a nematic or cholesteric phase for small D and usually stripe, triangular, and isotropic phases for large D. Furthermore, in the overviews some non-trivial phases that occur in between are shown. Details of the phases are discussed in the following subsections as denoted in the legend of figure 3.

Discussion of stable phases
In the following we discuss the stable phases that we observe in detail.
Note that in all snapshots depicting the phases the color denotes the density-like field ψ( x). The orientation field U ( x) is superimposed over ψ in the form of markers, which show the symmetry and direction. Furthermore, the size of the circle at the    markers denotes magnitude |U |. Note that the fields U and ψ; are determined with the same spacial resolution. However, we plot fewer orientation markers to keep the plots legible.
3.3.1. Isotropic phase For large D and large B l an isotropic phase is observed. In the isotropic phase the density is homogeneous and the magnitude |U | of the orientation field vanishes, i.e., there is no preferred orientation of the particles.

3.3.2.
Nematic and cholesteric phases For small D structures with strong orientational order occur. In case neighboring particles prefer to align i. e. for unmodulated interactions a nematic phase is found as depicted in figure 4(a). If neighboring particle prefer opposite orientations i. e. for modulated alignment interactions we observe a phase where a strong orientation changes continuously in a wave-like pattern. We call this structure a cholesteric phase.
While in an ideal nematic phase all particles possess the same orientation, in the cholesteric phase particles with similar orientation occur along stripes. Furthermore, the density is constant in the nematic phase but there is a stripe-like weak density modulation in the cholesteric phase.
In both phases, when defects are present, they tend to cluster in motifs that occur in close-by phases.  The orientation fields of the stripe phases depend only on n and are independent of (un)modulated alignment.

Stripe phases
For small B l and usually larger D stripe phases with a strongly modulated density occur. The stripe phases of modulated and unmodulated alignment are identical. For all odd n we observe the strongest modulations to appear where the density gradient is maximal as depicted in figure 5(a). In contrast, for all even n, figure 5(b), the orientational strength is maximal in the density maximum as well as in opposite direction (rotated by π/n) in the density minimum but vanishes at the flanks of the stripes.

Triangular phases
The triangular phases occur at larger B l than the stripe phases. Similar to the stripe phases, the orientation fields of the triangular phases depend whether n is even or odd; yet with exceptions when n matches a symmetry of the triangular lattice. Typically for even n the orientation is strong between neighboring density peaks, see Figure 6(b), while for odd n each density peak is surrounded by a ring of strong orientation as shown in Figure 6(c). In both cases the orientation vanishes at the density maxima, marking the triangular phases as plastic crystals. Exceptions of this general rule are n = 3 and n = 6: Figure 6  between triplets of neighboring density peaks. Figure 6(d) presents the triangular phase of n = 6, which is an oriented crystal. However, it only occurs as metastable phase as the dodecagonal quasicrystal is stable for these parameters. Figure 7 shows several rhombic phases for different n. The rhombic phases often appear at intermediate D between the nematic or cholesteric phase and other crystalline phases. In some cases the angles that occur might not be uniquely given: For modulated alignment next to the cholesteric phase, the rhombic phase often can be stretched or compressed at almost no energy cost, because there is a flat segment in the piecewise definition of f −k 2 in (14).

Rhombic phases
3.3.6. Honeycomb and Square phases The honeycomb and square phases, depicted in figure 8, add to the examples of plastic crystalline phases which appear in between the ubiquitous phases. Their topology is the same as reported for 2-fold particles in [15].

Dodecagonal quasicrystalline phase
We also find a stable quasicrystal, i. e. a structure with long-ranged order but without translational symmetry. The observed quasicrystal is shown in figure 9. It possesses 12-fold rotational symmetry and occurs for a system with modulated alignments and particles with 6-fold rotational symmetry. From a zoom-in to the quasicrystalline structure as shown in figure 9(b) one recognizes that along a ring around a local symmetry center there are 12 regions with strong orientation. The orientation in neighboring regions is rotated by π/6 as expected due to the modulation of the alignment interaction. Probably these regions are the features that stabilize the quasicrystals. As typical to quasicrystals, the Fourier transform of ψ exhibits not only the 12 main peaks but numerous quite strong satellite peaks as well. Therefore, the anisotropic interaction can indeed lead to the stabilization of lengthscales that differ from the length that is preferred by the free energy functional.
In principle a similar structure might occur for particles with 12-fold rotational symmetry and unmodulated alignment interactions. However, we have not found such a quasicrystal. Furthermore, we have not observed any stable quasicrystals with other rotational symmetry.

Defects and metastable states
Although we do not observe any stable quasicrystals with 5-fold rotational symmetry, it seems that 5-fold symmetry prominently occurs in defects along grain boundaries as shown in figure 10.
We also observe metastable quasicrystals like the one shown in figure 11 that are different from the dodecagonal structure reported in section 3.3.7.
The figures shown in figure 3 are not necessarily complete phase diagrams. In some cases it is hard to figure out whether a resulting structure is stable or metastable. For example in figure 3(b) some points in parameter space are marked with crosses. For these parameters we observe structures that seem to be mixtures of triangular, stripe, and possibly rhombic phase as shown in figure 12. Note that the corresponding pure phases possess a higher free energy than the mixture that we observe. Probably in these points the true equilibrium phase has not been found yet.

Conclusions
We have introduced a new PFC model for particles with n-fold rotational symmetry in two dimensions as it occurs, e. g., for patchy particles with symmetrically placed patches. Both the cases of attractive as well as repulsive patches have been considered leading either to alternating or the same orientation of neighboring clusters.
The PFC model is used to determine the phases that are stable for various rotational symmetries. We usually observe nematic or cholesteric phases in case a non-vanishing orientational order is preferred by the free energy. In the opposite limit we find stripe, triangular, or isotropic phases depending on how strong density modulations are supported. In between these phases, complex orderings with honeycomb, square, or rhombic symmetry occur. We even find a quasicrystalline phase with dodecagonal rotational symmetry. Therefore it is demonstrated that quasicrystals can be stabilized by interactions that only possess one length scale if in addition special binding angles are preferred. Note that we do not find the type of quasicrystals that has been reported to occur for patchy particles in simulations [25][26][27][28]. As we will discuss in the last paragraph this is probably due to lack of a hard core in our approach.
Since a lot of complex phases occur in our system, finding the global minimum might be hard and usually there are a lot of metastable states with interesting structures. At the end of section 3 we shortly comment on some examples including a metastable quasicrystal. However, the metastable states deserve more detailed  analyses in future works. For example, the domains that meet at grain boundaries prefer orientations that depend on the rotational symmetry and the orientational interactions. As a consequence the grain boundaries probably differ from the boundaries that are observed for isotropic particles. By adding noise the coarsening processes in systems with such grain boundaries can be explored, which we leave for future works.
The particles that we have in mind in this work do not possess a hard core. Therefore, the phases described here correspond to cluster crystals similar as in [10,[33][34][35][36][37]. In contrast, in many particulate systems a hard core prevents large overlaps of particles. Furthermore, a hard core can support the formation of some complex phases like quasicrystals [38][39][40]. Note that the quasicrystalline phases that have been observed in computer simulations of patchy colloids have been found in systems where the particles can hardly overlap [25][26][27][28]. As a consequence, in future we want to study mean field approaches with similar couplings between an orientational field and the density-like field as in this article. However, to model a hard core as well, the functional dependence on the density-like field has to be changed. A suitable way to describe the hard core is given by the so-called fundamental measure theory [41][42][43], which can be formulated in two dimensions [44] and which is known to lead to complex phases in case of the competition with an incommensurate substrate [45][46][47].

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.