Coarse-grained modeling of crystals by the amplitude expansion of the phase-field crystal model: an overview

Comprehensive investigations of crystalline systems often require methods bridging atomistic and continuum scales. In this context, coarse-grained mesoscale approaches are of particular interest as they allow the examination of large systems and time scales while retaining some microscopic details. The so-called phase-field crystal (PFC) model conveniently describes crystals at diffusive time scales through a continuous periodic field which varies on atomic scales and is related to the atomic number density. To go beyond the restrictive atomic length scales of the PFC model, a complex amplitude formulation was first developed by Goldenfeld et al (2005 Phys. Rev. E 72 020601). While focusing on length scales larger than the lattice parameter, this approach can describe crystalline defects, interfaces, and lattice deformations. It has been used to examine many phenomena including liquid/solid fronts, grain boundary energies, and strained films. This topical review focuses on this amplitude expansion of the PFC model and its developments. An overview of the derivation, connection to the continuum limit, representative applications, and extensions is presented. A few practical aspects, such as suitable numerical methods and examples, are illustrated as well. Finally, the capabilities and bounds of the model, current challenges, and future perspectives are addressed.


Introduction
The original phase-field crystal (PFC) model, introduced in 2002 [1], was developed as a simple way to incorporate elasticity and dislocations in continuum models in a manner similar to how interface and domain boundaries are introduced in traditional phase-field (PF) models. In the latter case, the predictions of PF models can be shown to be consistent in the asymptotic limit of vanishing interface widths with well-known sharp interface (SI) models [2] that explicitly track the position of a given interface subject to various boundary conditions (such as, e.g., the Gibbs-Thomson condition (GTC) for solidification or spinodal decomposition). PF models do not typically provide quantitative predictions on small length scales, i.e., on the scale of interfacial widths or suitable correlation lengths. Usually, their parameters are chosen to match the ones entering SI models [3][4][5] (e.g., the capillary length and coefficient of kinetic undercooling that enter the GTC). Similarly, PFC models do not quantitatively describe small length scale features, but in the appropriate limit they reduce to standard results. It is straightforward to show that in the long-wavelength limit, the PFC free energy reduces to traditional continuum elasticity theory [6] and that the dynamics incorporate vacancy diffusion [1,7]. It has been shown, numerically in two dimensions, that GBs can form spontaneously and their energy is consistent with the Read-Shockley equation [1,[7][8][9], that climb and glide of dislocations follow the Orowan equation [10], and in three dimensions that glide (climb) mediated sources of dislocation are consistent with Frank-Read (Bardeen-Herring) mechanisms [11]. More recently, it has been shown analytically that in PFC models the velocity of dislocations is determined by the Peach-Koehler force as expected in pure [12] and binary systems [13]. In addition, the predicted elastic fields around a dislocation agree quantitatively with continuum elasticity theory, encoding additional features such as anisotropies and non-linearities [14][15][16]. In many ways, the connection between PF and sharp interface approaches is analogous to the connection of PFC models with dislocation dynamics (DD) models [17][18][19], which explicitly move dislocation lines due to Peach-Koehler forces that are generated by the elastic field of other dislocations, defects, or externally applied forces. In particular, the coarsegrained PFC model referred to in the literature as amplitude expansion of the PFC, complex amplitude phase-field crystal or simply amplitude equations, on which this review focuses, allows a description of defects without resolving atomistic length scales, closely resembling the basic features of DD models. The advantage of this approach over DD is that dislocations and their main phenomenology appear naturally, following from the considered free energy functional. Therefore, no external rules would be in principle needed to determine the interaction, annihilation, or creation of any type of defect. At the same time, the method is not restricted to a single-crystal sample with pre-defined glide planes. However, it is worth noting that quantitative description of specific phenomena and materials would require an extended parametrization compared to minimal PFC-like models typically reported in the literature. Such extensions may be achieved with later formulations [20,21] but to date, they have not been explored extensively in this regard.
The complex amplitude phase-field crystal (APFC) model was originally derived by Goldenfeld et al [22,23] from the PFC model, which describes the evolution of the atomic number density during crystallization and the related dynamic processes [1,7,24]. While the PFC model can access diffusive time scales, the approach is limited by the need to incorporate density fluctuations on atomic length scales, thus requiring resolutions smaller than the lattice spacing. The main aspect of the APFC approach is to model the amplitude of the density fluctuations instead of the density itself. The idea of describing liquid/solid transitions by amplitudes that are real has been exploited in the past [25][26][27]. In Goldenfeld et al 's formulation [22,23], density fluctuations are described by complex amplitudes, η hkl , where hkl are Miller indices that describe specific crystallographic planes. The magnitude of η hkl is finite in a crystal and zero in the liquid state. Thus, it can be used to characterize a liquid-solid transition. Gradients in the phase of η hkl occur when the crystal state is strained, which provides information about the elastic energy stored in the crystal. In addition, the phase can describe the rotation of the crystal, allowing for the study of polycrystalline states (although, as noted in Sec. 5, there exist limitations). Finally, the combination of the magnitude and phase can describe dislocations in which large gradients in the phase do not lead to huge increases in the elastic energy as the magnitude of η hkl goes to zero. While the APFC model is formally derived from the PFC model, it is in principle possible to phenomenologically write down an APFC model as long as it has the correct long-wavelength behavior as has been done for PF models of various phenomena.
One of the most important features of the APFC model is that it provides a natural bridge between atomic and mesoscopic continuum length scales. In a single crystal state, the amplitudes vary slowly in space (depending on the orientation) but can be used to reconstruct the underlying atomic density fluctuations completely. On long length scales, it is straightforward to derive standard continuum elasticity through the phase of the amplitudes. Significant variations of amplitudes occur at defects and solid-liquid interfaces, still well describing the deformation induced in the lattice. The equations entering the APFC model, similarly to PFC, can be solved with simple numerical approaches. For example, using a uniform grid, Smirman et al [28] studied Moiré patterns in graphene films with the largest size system of 19.6 µm × 33.9 µm containing more than 25 billion unit cells (although it should be noted that these patterns contain no defects). When dislocations, grain boundaries, and interfaces appear, i.e. when a significant local variation of amplitudes occurs, more advanced numerical approaches can be considered to optimize the calculations. Indeed, these regions require the finest resolution, while a coarser one, typically much larger than the atomic spacing, can be used elsewhere. Adaptive meshing schemes then allow for simulation of large mesoscopic scales and at the same time completely retaining atomic information. Thus the APFC method allows simulations of atomistic features on continuum scales and should play an important role in understanding complex phenomena with multiscale features.
The rest of the review is organized as follows. Section 2 describes the original PFC model and the derivation of the APFC model. Section 3 outlines various numerical methods that have been developed to solve the APFC on regular and adaptive meshes. This is followed by Section 4 that provides a connection of the APFC model to traditional models of continuum elasticity and plasticity. Section 5 outlines the limitations of the approach and some extensions aimed at overcoming some of these constraints. Following this is Section 6 which describes some applications of the model to various physical phenomena. Finally, some conclusions and future outlooks are given in Section 7.
2. From phase-field crystal to the amplitude expansion 2.1. Origin of the phase-field crystal model The PFC model was proposed phenomenologically [1,7] to model elasticity and plasticity in crystal structures and can be written in terms of a dimensionless Helmholtz free energy functional, F , which is given as, and an equation of motion, ∂n ∂t = ∇ 2 δF n δn , where n is related to the atomic number density difference and ∆B 0 , B x , t and v are constants that may depend on temperature [24]. Although Eq. (1) can be derived [24,29,30] from the classical density functional theory of Ramakrishnan and Yussouf [31], the approximations used give rise to poor atomic-scale predictions in most materials since this free energy is minimized by an almost sinusoidal density fluctuations, while in metals for example n is very sharply peaked Gaussians at each lattice point. Nevertheless the periodic nature of the solutions of Eq. (1), which mimic a time average of microscopic atomic density [32] and evolves over diffusive time scales [33], make it useful for studying a large variety of physical systems such as multi-component polycrystals, liquid crystals, quasi-crystals and colloids as well as a broad class of phenomena including crystal growth and nucleation, heteroepitaxy, pattern formation, dislocation dynamics, grain boundary morphology and motion [7,[33][34][35][36]. PFC models have been developed also for less conventional materials and systems such as, for instance, active crystals [37][38][39][40][41], active colloids [42], and viral capsids [43]. The fact that the solutions are not sharply peaked means that they can be described by a few Fourier components. In this regard the density is written in terms of complex amplitudes, η hkl , as follows, where n o is the average density, G hkl = hq 1 + kq 2 + lq 3 are reciprocal lattice vectors, with q 1 = 2π(a 2 × a 3 )/(a 1 · (a 2 × a 3 )) and cyclic permutations of (1,2,3) the principal reciprocal-lattice vectors, and a j the vectors defining the primitive cell of the crystal lattice [44]. Note that the summation goes over both negative and positive G hkl 's with η −(hkl) = η * hkl such that n is a real field. In two dimensions, one may define G hkl as above with l = 0, q i = 2πRa j /(a i · (Ra j )) for i = j and R a 90 • rotational matrix (clockwise or anti-clockwise). All these definitions satisfy the condition a i · q j = 2πδ ij . Two illustrations of the quantities entering Eq. (3) in 1D are shown in Fig. 1, namely corresponding to a solid-liquid interface and a uniformly strained 1D crystal. Since PFC type models produce smooth solutions it is a good approximation to use the fewest number of complex amplitudes that are needed for any given crystal symmetry (see also Fig. 2). For example, only six η hkl (so three independent η hkl ) are needed for a 2D triangular lattice (more explicit examples are given in Sec. 2.3.2). G hkl entering approximations with the smallest number of modes are shown in Fig. 2. As discussed in the next section the goal of the APFC model is to derive equations of motion for the amplitudes.

Derivation
There are various methods for deriving the amplitude expansion from the original PFC model. Essentially, it requires a separation of length scales by assuming that the complex amplitudes vary on length scales much larger than the atomic spacing. In general this is the same assumption of all phase field models which require that interfaces or domain walls make a smooth transition from one phase to another. This is illustrated in Fig. 1 for a one dimensional liquid/solid interface for a system of atomic spacing a and interface width W . The "phase field limit" is such that a/W 1. For instance, for a twodimensional triangular lattice it can be shown [45] that in the limit that n o = 0 and the complex amplitudes are real and identical (i.e., η hkl = φ, for all hkl), they are described by traveling wave solutions (with velocity V ) of the form, where W is the width of the liquid/solid front which can be written [45] as where ∆B 0 ls = 8t 2 /135v is the value of ∆B 0 at liquid/solid coexistence and W m is the maximum value of W and is given by For ∆B 0 > 9/8∆B 0 ls no traveling wave solution exists as the solid is linearly unstable. Thus the phase field limit occurs when B x → ∞ and as such 1/B x can be used as a small parameter in a multi-scale calculation. In light of this, it is convenient to make the following rescaling, where τ = t/ √ vB x . Now the limit B x → ∞ corresponds to → 0. Goldenfeld and co-workers [22,23] report that to obtain rotationally invariant equations using multiple-scales analysis requires going to sixth order perturbations, which is an extremely tedious task, as to lowest order the resulting equations are not rotationally invariant. However, they have shown that this analysis gives the same result using a simpler renormalization group calculation. Other works addressed refinement and assessment of the general renormalization group approach [46,47].
To grasp the essence of the calculations without using these more rigorous methods, Athreya et al [23] developed a method that was coined "quick and dirty" that essentially obtains the same result in the W → ∞ limit. The basic idea is to assume that the amplitudes are constant on atomic length scales, i.e., u.c.
where u.c. is an integration over a unit cell and q is a sum over various G hkl . Since q is periodic in the unit cell, Eq. (8) is zero unless q = 0. This is a considerable simplification that reduces the number of terms that enter the free energy. For example, consider a term Only the first and last term for hkl = −(h k l ) give non-zero contributions using approximation Eq. (8), since they do not contain terms multiplied by a periodic function. Thus, in this approximation, Eq. (9) reduces to As discussed in the next section, contributions that arise from higher order polynomial terms will depend on the specific crystal symmetry under consideration. Terms containing the ∇ 2 operator are treated similarly noting that, assuming constant or slowly varying n o , Thus the Laplacian operator transforms as ∇ 2 → L hkl . While the effective operator on the right hand side of Eq. (11) appears to be anisotropic (due to the specific direction of the G hkl 's), it can be shown that the free energy is independent of the orientation of the pattern formed in n [48]. With these steps an energy functional which depends on amplitudes, F η , can be derived (see also Sec. 2.3). The dynamics of η hkl approximating (2) can be obtained by multiplying Eq. (2) by e −iGhkl·r and integrating over a unit cell, i.e., where V is the volume of a unit cell, which may be written as ‡ where the long-wavelength limit has been used in the last approximation. It is interesting to note that the equation of motion for the amplitudes are non-conserved, implying that an initial liquid (crystal) can completely transform in a crystal (liquid) locally. Nevertheless the density is a conserved quantity in a closed system and it is often important in liquid solid transitions since in liquid/solid coexistence the liquid and solid have different densities. In addition, the process of dislocation climb involves the mass (or vacancy) diffusion. In the original derivation of the APFC [22,23] the average density was assumed to be constant. The first inclusion of a spatially dependent density was reported by Yeon et al [49]. In this work n o was assumed to vary on the same length scales as the complex amplitudes and Eq. (3) should read Unfortunately, using the so-called "quick and dirty" method leads to an equation of motion for n o (and free energy) which contains terms like (1 + ∇ 2 ) 2 n and then implies that crystal state can be obtained from constant amplitudes or by a periodically varying n o (which of course violates the assumption the n o varies on the same length scales as the amplitudes). To overcome this difficulty several simpler models were proposed, which were shown to incorporate interfacial energy associated with the density difference at liquid/solid front as well as the well known Gibbs-Thomson effect [49]. The model can be written with dynamics The specific terms that emerge when averaged over a unit cell are discussed in the following section. This approach is also discussed in Huang et al [29]. If the amplitudes are assumed to be real (which eliminates the possibility of elastic and plastic phenomena) this reduces to Model C in the Hohenberg/Halperin [50] classification scheme that can be used to study phenomena such as directional solidification [51] or eutectic solidification [52,53]. Heinonen et al [54] use a similar free energy functional, but also incorporate momentum through the Navier Stokes equation and add the corresponding convective term to the dynamics of η hkl and n o . This has the advantage of including faster relaxation of elastic fields as discussed in Sec. 5.2. ‡ The functional derivative δF/δz * is computed treating z and z * as independent variables.

Formulas for amplitude equations
Let's consider the free energy Eq. (1) with constant average density n o and for the sake of simplicity the generic parameters The amplitude expansion is based on the approximation of n as from Eq. (3) with a finite set of M vectors G hkl , reproducing a specific crystal symmetry. This equation, exploiting that η −(hkl) = η * hkl , is here rewritten as where for simplicity G hkl is given a single subscript m and c.c. is the complex conjugate, highlighting the minimal set of amplitudes to be considered to approximate n. The free energy and the evolution law for the amplitudes can be obtained by exploiting the coarse-graining procedure introduced in Sec. 2.2, i.e. by integration over the unit cell of the phase-field crystal energy density (1), with n expressed through its amplitude expansion, Eq. (17) [48,[55][56][57][58].
To provide a general form of the free energy, consider separately the different powers of n entering Eq. (1), namely n k ({η m }, {η * m }) → ζ k . After averaging over a unit cell the following results emerge, with and neglecting terms including a factor K im+in with i = ±1, ±2 which would appear in ζ 2 and ζ 4 as G m with the same lengths are never parallel (or antiparallel), so K im+in = 0. Notice that terms as in the first sum in ζ 3 or the third sum of ζ 4 contributes if considering modes with two or three times the length of others, respectively (e.g. G 10 and G 20 in Fig. 2(b)). For a one-mode approximation of n through Eq. (17), i.e. by considering the shortest G m , and transformation (11), the excess term becomes u.c.
with G m = ∇ 2 + 2iG m · ∇ and L m = G m − |G m | 2 . In the one mode approximation, the length scales can always be re-parametrized such that |G m | = 1. Interestingly the term ζ 2 = Φ does not depend on the lattice symmetry, while ζ 4 can be written where ζ s 4 depends on lattice symmetry. Therefore, the free energy as function of amplitudes may be written The dynamics of the amplitudes, based on the PFC formulation in Eq. (2) and according to transformation (12) are given by where L m ≈ −|G m | 2 as in Eq. (13), and, from Eq. (18), 2.3.1. Multi-mode approximations. To model some crystal lattices, more than one mode is required in Eq. (17), i.e. more length scales are set through the choice of the reciprocal space vectors. In this case, ζ m reads as reported above, but the excess term takes different forms. However, it may be reduced to Eq. (20) through approximation [6,13]. For two lengths, R 1 = 2π/k eq 1 and R 2 = 2π/k eq 2 , corresponding to different lengths in the reciprocal space k eq 1 = 1 and k eq 2 = αk eq 1 , with α = 1 = k eq 2 /k eq 1 = R 1 /R 2 , the term including the differential operator in the dynamic would read [6] ( with and lengths have been scaled such that Therefore, the coefficient A can be rescaled by a factor α 4 /(α 2 −1) 2 and the same energy term as for the one mode approximation can be used. This result may be generalized for a lattice having N different length scales R = 2π/k eq and k eq /k eq 1 = α (noting k eq 1 = 1). Eq. (24) would read that for N l = 2, α 1 = 1 and α 2 = α reduces to Eq. (26). Then, under this approximation, Notice that in the presence of more than two modes, the coefficient of G 2 m cannot be taken outside the sum so it cannot be included in the coefficient A through rescaling as in Eq. (26).

2.3.2.
Results for specific lattice symmetries. Implementations of the APFC equations may be performed in a general fashion by considering Eqs. (18) and (23). This delivers a general framework suitable for changes in lattice symmetries and the number of modes used (eventually also different symmetries at once, see also Sec. 6.4). However, the specific equations corresponding to given lattice symmetries through the choice of reciprocal lattice vectors may be useful for analytic calculations and adhoc implementations. In the following, f s ≡ f s ({η m }, {η * m }) are reported for selected crystal symmetries used in literature, with the length of shortest reciprocal space vectors normalized to 1 (see, e.g., [6,[59][60][61] and Fig. 2). Triangular (TRI) symmetry (2D), one-mode approximation, N = 3: Triangular (TRI) symmetry (2D), two-mode approximation, N = 6: Square (SQ) symmetry (2D), two-mode approximation, N = 4: Square (SQ) symmetry (2D), three-mode approximation, N = 8: Body Centered Cubic (BCC) symmetry (3D), one-mode approximation, N = 6: Body Centered Cubic (BCC) symmetry (3D), two-mode approximation, N = 9 Face Centered Cubic (FCC) symmetry (3D), two-mode approximation, N = 7: Other symmetries may be considered, provided that the proper set of the reciprocal space vectors are known and that the encoded symmetry corresponds to a global energy minimum for some parameters (see Sec. 2.3.3). Alternatively, stability of phases/symmetries may be enforced with the APFC formulation outlined in Sec. 2.4.

Stability of phases.
In a relaxed, bulk crystal, real and constant amplitudes φ may be computed by energy minimization. For instance, for one-mode approximations and n o = 0, one gets the energy Letting ζ 3 = pφ 3 and ζ s 4 = qφ 4 where p and q where are integers, and minimizing the free energy given in Eq. (18), with respect to φ (δF [φ]/δφ = ∂h[φ]/∂φ = 0) gives the solutions, with ± the solution for C ≶ 0. For instance, for a triangular symmetry described by a one mode approximation (see Fig. 2 . Similarly, for a BCC lattice described by a one mode approximation (see Fig. 2) where M = 6, p = 48, q = 144 the result is . Moreover, the general stability of the solid phase described by a real amplitude φ 1,2 can be assessed by evaluating the condition F [φ 1,2 ] < F [0]. Notice that, F [0] is trivially 0 from Eq. (36), but it may have different values for n o = 0 as a non-zero average density would enter explicitly the energy (36) and modifies the value of the real amplitudes at equilibrium (see e.g. Ref. [6]). Phase diagrams can then be devised generally for both PFC and APFC approaches [6,62] by evaluating the relative stability of different phases described by φ. Generally, for a given set of parameters C and D, liquid phase results favored for values of B smaller than a critical value B c . This parameter phenomenologically encodes the role of the temperature. |B − B c | is often referred to as quenching depth. Notice that B c = 0 for C = 0.
When considering approximations with more modes, different values of φ should be considered for every set of amplitudes corresponding to different lengths of G m . Typically this task should be addressed numerically. Consider an approximation with K equal to the number of the modes of different length (under approximations introduced in Sec. 2.3.1). In this case the following function must be minimized, with M k the number of reciprocal space vectors for each considered mode (the solid arrows in Fig. 2). For instance, for the three-mode approximation of a cubic lattice in Fig. 2 are the symmetry-dependent polynomials resulting by substituting η j with the amplitude associated to the length of the reciprocal space vector they correspond to. To introduce an explicit example, consider the two mode approximation of the triangular symmetry (see Fig. 2(a)), i.e. {φ k } = [φ I , φ II ], M I = M II = 3, and ζ 3 (φ I , φ II ), ζ s 4 (φ I , φ II ) the polynomial resulting by setting η j = φ I for j = 1, 2, 3 and η j = φ II for j = 4, 5, 6 in Eq. (30). Plots of h(φ I , φ II ) for selected parameters (C = −2.0 and D = 1.0) are shown in Fig. 3. At a value B = 0.3 ( Fig. 3(a)), relatively close to the solid-liquid phase transition, the free energy has a single minimum corresponding to φ I ≈ 0.274 and φ II ≈ 0.087. By increasing the quenching depths, the global minimum shifts to φ I ≈ 0.215 and φ II ≈ 0.086 for B = 0.0. Moreover, another relative minimum appears (see Fig. 3(b)), which corresponds to a graphene-like phase. Some extended discussions on all the possible phases which can be described in two dimensions with combination of more modes can be found in Ref. [62].

Amplitude XPFC
A formulation based on the the so-called structural PFC (XPFC) [20,21], describing more detailed features and phenomena in crystalline systems such as, e.g. multicomponent systems, structural transformations, anisotropies, and extended defects [11,58,63], has been proposed in Ref. [58]. In a dimensionless form, the XPFC free energy F X reads where P and Q are parameters and X 2 (|r−r |) is the direct two-point correlation function at the reference density n o . In this approach, this function is typically expressed in the reciprocal space,X 2 (|k|). Following [58], it may be expressed as an envelope of Gaussian peaks associated with different modes of the periodic density or, in other words, to a family of planes of a crystal structure, [21] where w j controls the elastic and surface energies (the width of the j-th Gaussian peak), σ is an effective temperature parameter [64], p j and a j are the planar and atomic densities associated with the family of planes corresponding to the j-th mode, respectively, while k j is the inverse of the interplanar spacing for the j-th family of planes. Then, by assuming an amplitude formulation and volume average as in Sects 2.2-2.3, the polynomial in n that enters F X leads to terms similar to the energy in Eq. (15) except for the excess term which becomes [58] where the hat symbol denotes the Fourier transform, F −1 the inverse Fourier transform, andξ V an averaging (convolution) kernel in Fourier space that restricts the wave number to small values, approximately approaching the extension of the first Brillouin zone, which filters out spatial variations smaller than the lattice spacing. Interestingly, this model has been proposed with an ansatz for the amplitude expansion encoding different (two) lattice symmetries (see Sec. 6.4). This ansatz is expected to work with other forms of the energy and it consists just of a different formulation for Eq. (17) leading to results that may be formulated in terms of the equations reported in Sec. 2.3.

Numerical methods
In this section, two standard methods (finite difference and spectral) for solving first order in time partial differential equations that are applicable to APFC models are described. Following this, a finite element approach for solving APFC models is outlined and the description of a mesh refinement algorithm is reported.

Finite differences
In general there are many methods for solving an equations of the form where H(ψ) is a function of ψ. To solve it numerically it is useful to first consider integrating the equation over time from t to t + ∆t to obtain, The main question is how to approximate the integral in the above equation. In explicit methods only prior knowledge of ψ and its derivatives are used, i.e., where H(t) = H(ψ(t)). The simplest method, Euler's method, just retains the first term in the expansions, i.e., This approach must be supplemented by methods to evaluate spatial gradients in H, which in (A)PFC type models are typically even order derivatives, i.e., ∇ 2 , ∇ 4 , . . . . Often these are evaluated using a central difference formula. For instance, in two dimensions with a 5-points stencil (quincunx ), the Laplacian is given by where (x, y) = (i∆s, j∆s). Eq. (46), in conjunction with Eq. (45), is quite simple to implement for numerical integrations. Moreover, it is easy to incorporate different boundary conditions. However, the time step ∆t is limited by the grid spacing due to stability constraints, typically ∆t < α∆s −k , where k is the highest order spatial derivative (i.e., k = 6 for the PFC equation) and α is a constant that is model specific. If ∆t is too large, the solution very rapidly diverges (a pitchfork instability). The specifics of the origin of this instability are described in detail in Ref. [48]. It is possible to slightly reduce this instability by including next nearest neighbours as done by Oono and Puri [65]. This limitation is quite severe in PFC and APFC models as k = 6 in the former case and k = 4 in the latter. This instability can be avoided using semi-implicit approaches that are typically done in Fourier space. However, implicit or more generally semi-implicit approaches may be exploited, evaluating terms in the integrals in Eq. (44) within the range [t, t + ∆t], to have more stable numerical schemes (see also Sec. 3.3). Also, finite difference approaches may be combined with spatial adaptivity which may allow for efficient simulations (see also Sect. 3.4).

Fourier spectral method
Spectral methods solve differential equations treating variables as a sum of basis functions with coefficients to be computed, i.e., through a global representation. The so-called Fourier spectral method exploits the Fourier transform, typically in its discrete formulation for numerical integrations (therefore often referred to as pseudo-spectral, Fourier method). This method is particularly suited for periodic boundary conditions.
A key feature of this approach is that, in the Fourier space, differential operators become algebraic expression of the wave vector, e.g. ∇ 2 ψ(t) → −|k| 2 ψ k (t), where ψ k is the (discrete) Fourier transform of ψ. No finite difference approximations are then required if solving for ψ k (t), and ψ(t) may be then obtained through a (discrete) inverse Fourier transform. Moreover, efficient algorithms exist to compute ψ k from ψ and vice-versa, namely exploiting the Fast Fourier Transform (FFT) algorithm [71]. The adaptation of such approaches to phase field modeling in materials physics can be found in reference [72]. This method generally allows for splitting off the linear term in H and solving that part exactly, i.e., where L is a linear operator and N is a non-linear function of ψ. Indeed, in Fourier space, this would then read with N k the Fourier transform of N (ψ) and L k is an algebraic expression of the wave vector. Eq. (49) is an ordinary differential equation with solution Typically, the numerical instability in Euler's method occurs when L k is the most negative (i.e., at large wavevectors). However, in this method, e L k t is very small in this limit so that instability is completely avoided. To complete the picture, the non-linear term must be approximated as was done for H(ψ) in the preceding section. Considering Eq. (50) for ψ k (t + ∆t) and approximating (explicitly) N k (t ) ≈ N k (t) gives

Finite element method
The Finite Element Method (FEM) emerged as a particularly suitable framework for solving the APFC model's equations [16,60,77,78], besides being also employed in PFC studies in the first place [79][80][81][82][83]. Indeed, it conveniently discretizes partial differential equations (PDEs) while exploiting inhomogeneous and adaptive meshes. Within FEM, the PDEs are expressed in an integral form (weak form) over their domain of definition (Ω), typically having a rectangular/cubic shape. For the discretization of the resulting equations, a conforming triangulation T h of the domain Ω is considered, usually with simplex elements S ∈ T h (with characteristic size h).
In the context of APFC simulations, linear elements have been mostly adopted. This means considering a discrete function space of local polynomial of order 1 (P 1 ), namely h . To solve for complex functions, as η m , their real and imaginary part can be considered as two (real) independent unknowns. Alternatively, complex coefficients with real basis functions may be considered.
The FEM approach which has been used to solve APFC equations as in Eq. (22), features a splitting into two second-order equations for ∂η m /∂t and ρ m = G m η m (with m = 1, ..., M as in Sec. 2.3) [60,77]: This choice is convenient within the APFC framework as it allows the computing of relevant quantities straightforwardly as, e.g., the stress field, which may be rewritten in terms of both η m and ρ m and their spatial derivatives [16] (see also Sec. 4.2). Moreover, even though it is defined for G m , ρ m can be readily be used for computing L m , for instance when considering multi-mode approximations. From a numerical point of view, the splitting in Eq. (52) allows exploiting linear elements as only second-order operators appear, which translate to first order operators acting on elements of V 1 h in the weak form. With (f, g) := Ω f (r)g(r) dr the L 2 (Ω, R) scalar product, and considering the integral form of Eq. (52), the problem to solve then reads: ∀v ∈ V 1 h subject to an initial conditions η m (0) = η 0 m , and H({η})=∂f s /∂η m + Bη m + 3D(Φ − |η m | 2 )η m . The time derivatives are approximated by ∂a m /∂t = (a j+1 m − a j m )/∆t j and ∂b m /∂t = (b j+1 m − b j m )/∆t j , with ∆t j = t j+1 − t j the time step, and j ∈ N 0 the index labelling time steps. The time discretization is obtained through an implicitexplicit (IMEX) scheme. It consists of evaluating all the linear (nonlinear) terms in Eq. (53) implicitly (explicitly), i.e. at time t j+1 (t j ) [60,77], with a j+1 m , b j+1 m , c j+1 m , d j+1 m the unknowns to solve for. Eq. (53) consists of a set of nonlinear equations due to H({η}). This term can be generally linearized and handled through iterative approaches as Picard Iterations or the Newton method. A simple but effective approach, which can be exploited for methods introduced in previous sections too, consists of applying a one-iteration Newton method [60], i.e. approximating H(η j+1 ) as To solve Eq. (53), basis function expansions of unknowns are considered, e.g. a j+1 The FEM approach outlined above proved efficient in handling relatively large systems in both two and three dimensions, in combination with standard direct and iterative solvers within FEM toolboxes like, e.g., AMDiS [84,85]. Further improvements may be devised to increase the performances. An example is reported in [77] where the development of a dedicated preconditioner [86,87] allowing for fast solver convergence has been proposed and exploited for simulations of hundreds of nanometers domains in three dimensions for some materials.
The approach described in this section is also prone to coupling with other equations. Indeed, other variables would share spatial features with amplitudes. Coupling terms could be considered as additional terms entering ∂η m /∂t. At the same time, other equations may be discretized readily following the main FEM features described above (linear elements, operator splitting in second-order PDEs, IMEX time discretization). This has been exploited for instance when imposing mechanical equilibrium [16] (see Sec. 5.2), to simulate binary systems [13] (see Sec. 6.3), and to investigate the effect of magnetic field on small-angle grain boundaries [88].

Mesh adaptivity
Exploiting spatial adaptivity is a convenient strategy for performing efficient simulations with the APFC model [60,66,67,77]. Indeed, amplitudes are constant for relaxed crystals, oscillate with different periodicity according to the local distortion of the crystal with respect to the reference one (see, e.g., Fig. 1) and exhibit significant variation at defects and solid-liquid interfaces. Depending on the numerical approach and set of equations, one may devise different strategies to set a local refinement, e.g., based on error estimates or indicators.
An optimized local resolution based on the amplitudes oscillations, which works even for the standard approaches considered so far, has been achieved focusing on phases of the complex amplitudes, arg(η m ) = θ m . By looking at this quantity, it is possible to determine the wavelength of oscillating amplitudes λ m [77]. Then for a good resolution of all the amplitudes, the discretization h should be a fraction of the smallest λ m , i.e. h amp = min m (λ m )/n, with n ≥ 10.
To use this criterion in practice, the deformation, strain and/or rotation fields must be derived from amplitudes. This will be discussed in detail in the following section (see Sec. 4.2). In addition to the oscillation of amplitudes, a refinement for the interfaces and defects controlled by h min where |∇Φ| is significantly larger than a relatively small threshold ς and imposed as finest resolution in the mesh is considered [60], while a large discretization bound h max is defined for region where Φ ∼ 0 or where θ m → 0 (i.e. for constant amplitudes). Summarizing these concepts, this method ensures a local discretization, h, as This approach has been exploited together with the FEM approach outlined in Sec. 3.3, in particular within the FEM toolbox AMDiS [84,85]. However, it is expected to work with any real-space method readily. Further optimization of the mesh refinement can be achieved by a polar representation [66,67] which involves, however, some changes in the amplitude equations, the coupling with additional fields, and other technical details to be considered. An examples of an APFC simulation performed with the adaptive refinement strategy here outlined is given in Fig. 4.

Elasticity
The elastic properties in the amplitude expansion arise from the term A M m Γ m |G m η m | 2 (see Eq. (28)). Indeed, all the other terms in the free energy do not give rise to gradients in the phase of the amplitudes and as such do not contribute to the elastic energy. To obtain the consequences of this term it is useful to consider deformations (u ≡ u(r)) from a perfect lattice, i.e., where θ m ≡ G m · u and φ m is weakly dependent on u (see a 1D illustration in Fig. 1(b)). This leads to where in the last line higher order gradients in u have been neglected. So that where u ij ≡ ∂u i /∂x j , G m i is the i-th component of G m and the Einstein summation convention is used. Eq. (58) contains linear and non-linear terms. In terms of the non-linear Eulerian-Almanasi strain measure (U) [57,73] with elements §, § The strain measure U belongs to the general class of strain (material, Lagrangian) called Seth-Hill tensors ε n = (1/n)(C n − 1), with C = F T F, F ij = ∂x i /∂X j the deformation gradient and x and X the spatial (eulerian) and material (lagrangian) coordinates respectively, such that dx = FdX and dX = F −1 dx [89][90][91][92][93]. U corresponds to ε −1 . This definition mixes a Lagrangian tensor due to the dependence on F T F (an Eulerian tensor would depend on FF T ), with an Eulerian strain measure 1 − F −1 (a Lagrangian strain measure would depend on F − 1), see also Ref. [73].

Eq. (58) can be written as
The elastic part of the free energy is then The components of the stress tensor defined as where C ijkl is the elastic modulus tensor [94] are then given by Thus Eq. (63) provides a general formula for the elastic moduli for arbitrary crystal symmetry. Some specific examples are given below.
Examples: For a free energy with a single mode, i.e., containing the term n(1+∇ 2 ) 2 n/2, 2D triangular and 3D BCC structures minimize the free energy in certain parameter ranges. At a minimum these systems can be described by modes with the same length scale and thus Γ m = 1 and φ m = φ, ∀m. Following the definition of G m as in Sec. 2.3.2 for these symmetries (one-mode approximation), Eq. (61) gives For the FCC symmetry in the two-mode approximations (see Sec. 2.3.2), Γ m = 1/16 , ∀m. This gives where η m = φe −iθm for i = 1, ..., 4 and η m = ψe −iθm for i = 5, ..., 7.
One of the difficulties in parameterizing PFC models is that the ratio of the elastic moduli cannot be changed in the one mode triangular and BCC cases. However, it is interesting to note that in the FCC case, the ratio of the elastic moduli depends on ψ, which in principle can be tuned. It suggests that adding more length scales will allow for more tuneability in the models as shown in XPFC models [21]. However, it is important to note that if the added vectors have the same symmetry as the original ones this will not change the ratios.

Strain and stress field from the amplitudes
When examining the results of APFC simulations, it is useful to develop methods to extract the strain and stress fields directly from the complex amplitudes. As shown by Salvalaglio et al [14] the displacement field, u that enters continuum elasticity field can be extracted directly from the phase of the amplitudes (θ m ). In two dimensions (2D), inverting Eq. (56), the expression is with (i, j) = (x, y) and cyclic permutations, ij is the 2D Levi-Civita symbol, l and m label two different amplitudes,p =x ×ŷ the normal vector of the xy-plane and θ m = arg(η m ) = arctan [Im(η m )/Re(η m )]. In three dimensions (3D) it can be shown that with (i, j, k) = (x, y, z) and cyclic permutations, and l, m, n, labelling three different amplitudes. These quantities are discontinuous. However the component of the (linear) strain tensor U L become expressions of ∂θ m /∂x i with which is continuous almost everywhere in the solid phase, with a singularity for vanishing amplitudes in correspondence of phase singularities, e.g., at the cores of defects. Then, with a regularization for these amplitudes (see also Sec. 4.4), elastic field can be readily computed and conveniently exploited. In two dimensions, for U L and the rotation field ω = ∇ × u we then get Explicit expressions for 3D strain and rotation fields can be found in Ref. [14]. The stress field can then be computed through the Hooke's law (62). In 2018 Skaugen, Angheluta and Viñals [12] derived an expression for the stress tensor, σ ij from the density field using the standard definition of σ ij , i.e, where ∆F = F (n(r + u)) − F (n(r)) and u is the displacement field. This gives where P = f − n(δF/δn) is a pressure term summing up to the mechanical stress, with f the integrand in Eq. (1), the second term arising when considering mass-conserving deformations [95], and L ≡ 1 + ∇ 2 . In terms of amplitudes, integrating over the a unit cell with n expressed via Eq. (17) and neglecting the pressure terms gives [16] for one-mode approximations, while it can be generalized for more modes accounting for the full L m operators (see Eq. (20)).

Plasticity and defect dynamics
As seen in previous sections, the amplitude formalism can describe the elastic behavior of crystals as encoded in the PFC model. Moreover, by focusing on singularities in the corresponding phases, the motion of defects may be connected to the evolution amplitudes [12,13,15,96]. A dislocation in a crystalline lattice corresponds to a discontinuity in the phase θ m . At the same time, a dislocation with Burgers vector b is defined by du = b [97], thus it can be shown that dθ m = −G m · b = −2πs m , where s m is the winding number. As discussed in Ref. [12], a vortex solution for amplitudes at dislocation cores may be assumed, that reads η m ∝ x − is m y with s m = ±1. The Burgers vector distribution of a dislocation can be defined as a localized (vectorial) topological charge bδ(r − r 0 ) with r 0 the nominal position of the dislocation core, assumed pointwise from a continuous point of view. By extension, the Burgers vector density can be defined to be , with d indexing the dislocations and D their total number. To connect this quantity to amplitudes, note that the position of the core is where the amplitudes go to zero. Therefore, following the theoretical framework reported in [98][99][100], a change of coordinates from the canonical one to the amplitudes' components can be considered. Namely, for point dislocations in two dimensions, or straight dislocations in three dimensions, one gets with D m the Jacobian determinant of the coordinates' transformation, [12,98,99], and implying the Einstein summation convention. Aiming at the velocity of dislocations, the dynamics of B(r) is considered. Exploiting that the determinant fields D m have conserved currents [100], and that a similar continuity equation holds true for δ(η m ), from Eq. (73) the equation of motion for B i may be written, where the last term was obtained by transforming back the delta function to spatial coordinates. For dislocations moving at a velocity v d , it also follows that . Therefore, by equating this latter expression with the corresponding quantity in Eq. (75), the dislocation velocity can be related to the evolution of amplitudes as At the dislocation core, a few simplifications may be considered. For the amplitudes which are zero at the dislocation core, while others do not contribute to Eq. (76). The latter term in Eq. (77) is obtained by imposing again a form for amplitudes as in Eq. (56) and retaining the lowest order only in φ m and θ m . Combing all the equations reported above gives where U ij = (∂ i u j +∂ j u i )/2. This equation is consistent with the classical Peach-Koehler force [97]. For the case of a 2D triangular lattice or a 3D BCC crystal where it is possible to construct the lattice by retaining only one mode of the lowest order (with |G m | = 1, Γ m = 1), the velocity takes the form with M a mobility factor. With this formalism, the dynamic of defects may be obtained once ∂η m /∂t are known. This applies independently to the specific contributions affecting the dynamics of amplitudes. See, for instance, an application to binary systems in Sec. 6.3. The equations presented here apply for point dislocations in two dimensions or straight dislocations in three dimensions. A generalization to curved dislocations in three dimensions has been recently introduced in Ref. [96].

Comparisons with elasticity theories
As noted in previous sections, the APFC model may be employed to the study elasticity and plasticity in crystalline systems. A few prototypical cases have been investigated, delivering direct comparisons with predictions from other theories [14,16,101]. Of particular note is the comparison with continuum elasticity results, as the coarse-grained nature of APFC may deliver advanced/improved continuum approaches.
A representative case is the elastic field generated by dislocations at mechanical equilibrium, which is well known in the continuum (linear) elasticity for isotropic media [97,102]. In the APFC model, configurations with dislocations in prescribed positions may be obtained with different approaches. The phase of amplitudes σ hkl can be initialized with singularities as discussed in Sec. 4.3 at given positions and then the APFC model is used to minimize the free energy. By restricting the description to 2D crystals for the sake of simplicity, a convenient approach consists of setting phases the displacement field of an edge dislocation having Burgers vector b = bx and ν the Poisson's ratio [97]. Alternatively, an initial strain that induces the formation of dislocations can be considered. For instance, a pair of dislocations having the Burgers vector ±b edge is obtained by defining layers with initial deformation u = [Dx, 0] with D = ±b/L and allowing the system to relax [60]. Dislocations move when Peach-Koehler force is finite assuming no barriers exist (see Sec. 5.4). As discussed in Sec. 5.2, for dynamical configurations, corrections are needed to account for mechanical equilibrium within (A)PFC. Special cases are the configurations where defects do not move, and relaxation given by dynamical equations effectively approaches mechanical equilibrium. These may be represented, for example, by equally spaced arrays of dislocations alonĝ x andŷ with alternating Burgers vectors, i.e., a "grid" where four defects with the same Burgers vectors surround another one with opposite Burgers vector. It is worth mentioning that a single dislocation, in the absence of external stress, would be in principle stationery too (as the Peach-Koehler force is zero). Still, its elastic field would inherently interact with the boundaries of any finite simulation domain as it is longrange, with energy dependent on the system size and diverging for an infinite medium. A possible solution would be studying a single dislocation in a finite crystal [16], which, however, is expected to induce changes in the elastic field [97,103,104]. Fig. 5 shows the elastic field of a dislocation belonging to a two dimensional grid with alternating Burgers vector alongx andŷ. Both strain components resulting from computing Eqs. (69) (Fig. 5(a)) and stress components from Eq. (72) (Fig. 5(b)) are shown. These fields agree well with the field expected in classical continuum elasticity [97]. The elastic field obtained from Eqs. (69) is to some extent easier to compute as it involves only the first derivatives of amplitudes. Still, they are singular at the core of vanishing amplitudes, here regularized by setting to 1/(|η m | 2 + δ) with a small δ as prefactor in Eq. (68). On the other hand, the elastic field from Eq. (72) does not require such a numerical regularization. This approach involves higher-order derivatives than Eq. (68), which can be handled efficiently when combined with a proper splitting of the APFC equations (see also Sec. 3).
More insights are given in Fig. 5(c) and Fig. 5(d). Therein, a comparison of the stress field components obtained with different continuum theories for representative isolines (panel c) and along lines crossing the defect core (panel d) is reported. In particular, it shows the stress fields components computed from the APFC simulation, namely Eq. (72) and Eq. (62) with U L ij from Eq. (69) with φ 2 = 3 m=1 |η m | 2 /3. These fields are compared with the non-singular isotropic theory (NS) reported by Wei Cai et al in Ref. [102], and σ NS zz = ν(σ NS xx + σ NS yy ), with σ 0 = Eb x /(4π(1 − ν) 2 ), E the Young modulus, ν the Poisson ratio, and c a parameter controlling the extension of the core-regularization (c = 0 reduces to classical continuum elasticity (CE) formulations σ CE [97]). The triangular symmetry considered here, which results isotropic, and under the plane strain condition, gives µ = λ = 3φ 2 while E = µ(3λ + 2µ)/(λ + µ) = (5/2)φ 2 , and ν = λ/(2λ + 2µ) = 1/4 ¶. Another comparison with continuum elasticity is provided with a regularized formulation of the stress emerging in the framework of strain-gradient elasticity (Helmholtz type) [105,106] σ GE and σ GE zz = ν(σ GE xx + σ GE yy ), with K n (r/ ) the modified Bessel function of the second type, and a characteristic internal length parameter of the material. The elastic field obtained from APFC simulations encodes a smoothing similar to the non-singular theories in Eq. (81) and Eq. (82). A good agreement is obtained with c = 2a 0 and = a 0 . However, notice that these parameters are expected to vary for different quench depths as they are related to the extension of the core [102,105] and this shrinks with decreasing the temperature. It is worth mentioning that strain gradient terms may be indeed identified in Eq. (57), supporting the qualitative agreement shown in Fig. 5. For isotropic materials, a more accurate description is actually given by the so-called Mindlin's isotropic first gradient elasticity, which feature two characteristic lengths [107][108][109] and may therefore provide descriptions closer to the APFC results. Comparisons for 3D configurations and for rotation fields from Eq. (69) can be found in Ref. [14].

Large tilts: the problem of beats
Complex amplitudes consistently describe deformations, i.e., the energy is rotationally invariant while accounting for elastic energy associated with distortion with respect to the reference state (see Sec. 4.1). However, the larger the rotation with respect to the reference crystal (described by Eq. where G m i (Θ) = G m j R ij (Θ) and R ij (Θ) is the counter-clockwise rotational matrix. Therefore, oscillations of η Θ m have a wavelength 2π/|∆G m (Θ)|. This leads to a crucial two-fold limitation for the APFC model. On one side, the spatial resolution required to discretize the corresponding equations depends on their relative orientation with respect to the reference lattice encoded in G m . For large rotations this results in significant variations of the amplitudes over lengths approaching the lattice spacing, inconsistent with the assumption in their derivation and also requiring mesh sizes approaching the ones required in the PFC model. On the other side, while the energy of a single crystal remains rotationally invariant, the rotational symmetry of bicrystals is lost, and unphysical grain boundaries are obtained for large relative tilts corresponding to small or no deviations in the density field n (e.g., when rotating a 2D triangular lattice by ∼ 60 • ). An illustration of this behavior is reported in Fig. 6. When increasing the relative rotation of a circular inclusion, the oscillation of amplitudes increases requiring finer mesh as illustrated by Re(η 1 ). Even though the fields are properly resolved, unphysical grain boundaries appear in Φ for θ 30 • (e.g., according to symmetry, θ = −10 • and θ = 50 • should coincide, as well as θ = 60 • should have no defects with a Φ uniform).
An attempt to overcome this issue followed the first publications on the APFC model and consists of a polar representation of amplitudes [66]. In practice, the complex amplitudes are expressed in terms of the real fields φ m = |η m | and θ m = arg(η m ). The resulting set of equations for ∂φ m /∂t and ∂θ m /∂t derived from Eq. (22), have issues related to the discontinuous nature of θ m and that φ m vanishes in the liquid phase, in principle requiring robust and structured regularization algorithm. Therefore, further approximations are introduced [66]: i) a hybrid formulation exploiting the aforementioned polar representation only for crystal bulk, i.e. away from defects and interfaces, while solving the equations for the complex amplitudes everywhere else; ii) neglecting third and higher-order spatial derivatives of φ m and θ m in their dynamics and iii) assuming that gradients in the phase are zero within grains. This method has been shown to allow for efficient inhomogenous spatial discretization for numerical methods working in real space.
Recently the same issue has been addressed by exploiting a Cartesian representation of the amplitudes and allowing for local rotation of the basis vector G m [67,68]. This model considers a set of locally rotated amplitudesη m such as η m =η m e −i∆Gm(Θ)·r . A rotation field Θ is then computed such that η m have vanishing oscillation, i.e., satisfying the condition The local rotation field may be explicitly extracted from amplitudes, e.g. exploiting the results reported in [14]. Then, it may be shown [67,68] that operators defined in the rotated system, O Θ , applied to rotated fields, f Θ , transform as The evolution for η Θ is evaluated while computing G m (Θ) everywhere. This approach still requires a proper numerical implementation [67], but has been proved successful in describing crystal structures through the "rotated" amplitudes avoiding beats due to crystal rotation, exploiting efficient mesh refinement (see Sec. 3.4), and matching the dynamics obtained by the original amplitude expansion. Importantly, this approach has also been combined with an algorithm selecting the closest reference crystal for a given local orientation [68] which avoids the presence of unphysical grain boundaries, at least in two dimensions for triangular lattices.

Elastic relaxation and mechanical equilibrium
The dynamics of the PFC model and, in turn, its amplitude expansion approximation, was initially assumed to be overdamped, i.e. driven by minimization of the corresponding free-energy functional through a gradient flow [1,7]. Although this setting can be justified in some circumstances, it constrains the dynamic to diffusive timescales. This may lead to some issues for the description of elastic relaxation, which usually occurs on faster timescales with respect to the diffusive dynamics of the density field. A few investigations addressed these issues, delivering either a framework able to ensure mechanical equilibrium at every time, describing the limit of instantaneous elastic relaxation [15,16,75], or modeling explicitly elastic excitations [54].
In the work of Heinonen et al [75,115], the amplitudes are expressed similarly to Eq. (56), assuming small displacements in u. Then a formal separation of the timescales of the field φ m from the field θ m , is considered. To ensure mechanical equilibrium, i.e. ∇ · σ = 0, it is then demonstrated to be equivalent to solving at every step after solving for ∂η m /∂t. In [75], a factor φ −2 m appears in the second-last term in (87). However, as discussed in [115], this expression allows for a more formal connection to the displacement u. Moreover, equilibrating Eq. (87) would corresponds to a real energy minimization problem.
A different approach, which computes the mechanical equilibrium deformation from the incompatible one, fully accounting for the singular distortion of defects as conveyed by n and/or η m has been proposed in Ref. [15] for PFC and then translated to APFC in Ref. [16]. Therein, the smooth distortion u δ i required to fulfill mechanical equilibrium is determined, and then the amplitudes are corrected as η m.e. m = η m e −iGm·u δ . In brief, the smooth stress, σ δ ij , to be added to the stress field computed from the amplitudes, σ η ij (see also Sec. 4.2), to satisfy mechanical equilibrium is obtained through the Airy Function (χ) formalism: where B(r) the Burgers vector density, and ν, λ and µ as in Sec. 4.4, while u δ is then computed exploiting a Helmholtz decomposition into curl-and divergence-free parts, Once u δ i is calculated, correction to the amplitudes can be imposed. This approach has been shown to work well in two dimensions for isotropic materials, while its generalization to three dimensions is non-trivial due to the Airy function formalism. A more general method to correct n by computing u δ in three dimensions has been recently proposed in Ref [96] for PFC, and it is expected to work for the APFC model.
In Ref. [54], a model accounting explicitly for elastic relaxation has been considered by coupling the mesoscale description of the microscopic structure of the materials achieved by amplitudes to a hydrodynamic velocity field. It recovers the instantaneous relaxation as a limit of the model. It consists of describing the crystal lattice through η m and a slowly varying density field, n o , via the energy (15). The evolution laws are then derived accounting for mass density and momentum density conservation and read with v the velocity field, Dv/Dt = ∂v/∂t + v · ∇v, Q m = ∇ + iG m , and µ η , µ n , µ B , µ S are parameters. Previous attempts to include fast time scales in the dynamics introduced an explicit second order time derivative in the equation of motion for the PFC mass density field [116,117]. This approach gives rise to short wavelength oscillations accelerating relaxation processes, but fails to describe large scale vibrations [55]. The model described by Eq. (90) gives the correct long wavelength elastic wave dispersion relationship (ω ∼ k).
A key test case for all the approaches reported in this section is the shrinkage of rotated grains (see Fig. 7). Their results consistently show a faster dynamic in the limit of instantaneous mechanical equilibrium [12,16,75] while tuning of parameters in the model reported in Eq. (90) allows for the investigation of intermediate regimes [54].

Control of interface and defect energy
The original APFC (or PFC) model contains a small set of parameters which limits quantitative fitting to match experimental measures or theoretical calculations. In Ref. [60], it has been shown that the addition of a single term to the free energy functional can be used to control the solid-liquid interface and defect energies in a well-controlled fashion, without affecting the crystal structure. Exploiting the information conveyed by Φ = 2 M m |η m | 2 , which is a measure of the crystalline order, and in analogy with the gradient term of order parameters in interfacial free energies [119], an additional energy contribution can be phenomenologically introduced in Eq. (21), reading where β is a free parameter. This leads to an additional term to Eq. (22) as For small β, this additional contribution is found to change the interface and defect energy linearly with β, while deviations are observed for large values. Fig. 7(b) shows the tuning of symmetric tilt grain boundary energies by β due to the local change in the defect-core energies [60]. Notice that, due to the issues discussed in Sec. 5.1, it is not possible to compute the whole range of θ only by increasing the relative angle (see also [9]). In this case, energy values for theta ≷ 30 • are obtained with two different simulation settings. The framework reported in [68] would allow addressing these calculations without considering such different settings.
It is worth mentioning that formulations allowing for tunable energies at defects and interfaces similar to the one discussed here can be devised from microscopic length scales exploiting smoothing kernels in Fourier space [120,121].

Lack of barriers
In the derivation of the amplitude equations it was implicitly assumed that the atomicand meso-scales (interface widths, etc.) completely decouple. It appears that this approximation eliminates barriers for defect or grain boundary motion. Huang has shown that incorporating the first-order coupling of the atomic and mesoscales leads to interface pinning [118]. Consider multiplying the equation of motion by e −iq·r and integrating over a unit cell while keeping terms previously assumed to be zero. This leads to additional terms in Eq. (22). For instance, for a triangular lattice: where A u.c. is the area of a unit cell and with (· · · ) implying six other similar terms that contain a e −iq·r term (see reference [118] for details). The last term(s) in Eq. (93) implicitly couple atomic (e −iq 0 y ) and slow scales (η m ) terms. The equation for the average density becomes To understand the consequences of this coupling, Huang derived an equation of motion for a liquid/solid front moving in the y direction with slow variations in the x direction using the projection operator method of Elder et al [5]. In this method a coordinate transformation from (x, y) to (u, s) is made where u is a coordinate normal to the interface position and s is parallel. Equation (93) (in the limit L m ≈ −|G m | 2 = −1) is multiplied by ∂η m /∂u and Eq. (95) by ∂n o /∂u and integrated over u in the inner region. In the outer regime the Equations (93) and (95) are linearized around a liquid state and then solved using Green's functions. The inner and outer solutions are then matched such that the chemical potential is continuous across the interface.
One main result of these calculations is the equation for the interface normal velocity, v n , given by where c 0 is the kinetic coefficient, λ ∝ ∆n 0 0 δµ(0, s), ∆n 0 0 is the difference in liquid/solid density, δµ(0, s) is the chemical potential difference from equilibrium along the interface, γ is the surface tension, κ is the curvature, p 0 is the pinning strength, h is the distance from the front and φ is the phase. Expressions for each of these terms is given in Huang [118]. This equation coupled with mass diffusion in the outer regions (η m at equilibrium liquid values) and the usually matching condition v n ∆n 0 0 = ∂δµ/∂u| 0 − − ∂δµ/∂u| 0 + constitutes a free boundary problem.
If gradients in h are assumed to be small, Eq. (96) reduces to In the limit of non-conserved dynamics (fixed λ) this is a driven sine-Gordon equation introduced by Hwa et al [122] to study, when thermal fluctuations are included, the interface roughening during crystal growth. Huang showed that the pinning term can lead to step by step growth of the interface as is observed in experiments and even completely arresting the growth if the driving force (λ) is too small, as illustrated in Fig. 7(c). It is also shown that the pinning strength increases as temperature (controlled by B = ∆B 0 ) or the elastic moduli (controlled by A = B x ) are lowered as both have the effect of decreasing the width of the liquid/solid domain wall. Later, Huang [123] extended this work to a binary system with a eutectic phase diagram and derived more general expressions for the surface energy and barrier strength as a function of concentration, temperature, and crystallographic orientation of the liquid/solid front.

Solid-liquid interfaces and the phase field limit
Solid-liquid interfaces are regions where n may vary over length scales larger than the atomic spacing. Therefore, the APFC model may be exploited to focus on these regions while neglecting the fine details at the atomic scale elsewhere [124]. Real amplitudes have been first considered to address the modeling of solid-liquid interfaces in the seminal works by Khachaturyan [25,26]. Therein, the order parameters resemble the ones entering classical phase-field approaches [48,[125][126][127] and they may be linked to atomistic descriptions. They can be used, for instance, to account for bridging-scale descriptions of elasticity effects by means of additional contributions as, e.g., in the presence of precipitates, alloys, or point defects. [128][129][130][131][132]. However, this approach does not directly encode rotational invariance and elasticity associated with the deformations of the crystal lattice. In Refs. [45,61,133], traveling waves characterized by the ansatz (4) have been shown to describe the solid-liquid interfaces within PFC quite well near melting. Real amplitudes result in a classical phase-field model. Indeed, it is shown that a general form for the free energy can be obtained by considering real amplitudes, where the parameters a, b, c, d depend on the lattice symmetry and the number of modes considered. Different crystalline cubic lattices, and their effect on growth dynamics are still retained [61]. In addition, the framework is consistent with atomistic simulations and can be used for matching parameters to specific materials. In Refs. [124,134] similar underlying ideas led to a phase-field model connecting anisotropic surface energy and corresponding Wulff shapes to the lattice symmetry of various crystals through the choice of reciprocal lattice vectors. The model remarkably encodes a regularization term leading to corner rounding of faceted shapes similarly to diffuse interface theories [135][136][137]. Amplitudes are assumed to be real, but they are still considered separate variables. In the notation adopted in this review from Eq. (21), and assuming zero average density, this gives  99) is similar to Ginzburg-Landau free energies entering multiorder-parameter phase-field models. The higher-order gradient contribution [∇ 2 φ m ] 2 enforces the rounding of corners appearing among facets. A coefficient may be also introduced to tune its influence [134].

Grain growth with dislocation networks and small-angle grain boundaries
The PFC model has been exploited to investigate rather small systems due to the atomic-scale resolution. According to the features described in Sec. 4 and 5, the APFC is especially suited to describe systems with small deformation and rotation while including isolated defects such as dislocations. Examples include small-angle GBs in graphene structures [9], GBs premelting and shearing in BCC iron [139], and the dynamics of small-angle GBs in general [73]. In two dimensions, it is possible to examine systems on the micrometer scale [28,77] (see, e.g., Fig. 8(a)). A recent, remarkable application at this length scale is the simulation of sub-boundaries formation due to orientational gradients in thin aluminium films [76,140] (Fig. 8(b)).
The limitation in size for PFC becomes even more evident in three dimensions, requiring advanced numerical methods to simulate rather small systems [10,87]. The APFC model has been proved powerful in addressing the study of defects in crystalline systems in three dimensions [14,77,138]. In particular, small-angle grain boundaries can be well captured and also characterized thanks to the advanced description of elasticity as described in Sec. 4. Representative cases are the shrinkage of dislocation networks forming at the boundaries between rotated inclusions and unrotated surrounding matrix (see Fig. 8(c)), also in combination with additional effects (see also Sec. 6.3), and the growth of slightly misoriented crystal seeds (see Fig. 8(d)). Interestingly, the shrinkage or rotated inclusions and the resulting dislocation networks have been proposed directly using a classical PFC approach [10]. This investigation delivered very similar results to the ones obtained by APFC, as reported for instance in Fig. 8(c), thus assessing the coarse-graining achieved by the APFC model in an applied case. The shrinkage of grains is generally associated with their rotation. A fingerprint of this process emerges in APFC, as shown in Ref. [14] where rotations are tracked thanks to Eq. (69). Therein it is shown that when defects at the boundary of a grain get closer, their deformation fields superpose, increasing the effective orientation of the grain.

Binary systems
Coarse-grained approaches are often required in multiphase systems and alloys to handle simultaneously the deformation induced in the lattice, the resulting phase separations leading to Cottrell atmospheres [141][142][143], and effects on dislocation motion. The APFC model has been proved powerful in describing these effects at the mesoscale for binary systems, beyond results achieved by focusing on either atomistic or continuum length scales [144][145][146][147][148][149]. Also, it can be used to study these systems comprehensively, without focusing on concentration profiles, stress distribution around dislocations, and the forcevelocity curves for defect motion separately.
The original binary PFC model [24] is formulated in terms of the dimensionless atomic number density variation field and a solute concentration field ψ. In the APFC model, the expansion Eq. (17) is considered and a Vegard's law for the lattice spacing R = R 0 (1 + αψ) is assumed with α the solute expansion coefficient. This results in an energy [6,13] with definitions as in previous sections and w, u, Y, K, are additional model parameters as described in Ref. [24]. Dynamics in terms of ∂η m /∂t is then described by Eq. (13) with energy (100) and ∂ψ/∂t = ∇ 2 δF αψ /δψ, similarly to (16). It can be shown that, given G m the basic wave vectors corresponding to a pure system, the equilibrium wave vectors for binary systems read G eq m = G m √ 1 − 2αψ [29]. This approach allows the study of solute segregation and migration at grain boundaries, eutectic solidification, and quantum dot formation on nanomembranes [6,13,74,150]. A similar approach has been exploited to accurately describe the interactions among grain boundaries and precipitates in two-phase solids [59,69].
By applying the framework illustrated in Sec. 4.3 to this model, the velocity of dislocations including effects of the solute segregation has been also derived. By retaining only one mode of the lowest order (with |G m | = 1) and using the expression for ∂η m /∂t for binary systems into Eqs. (74)- (76) Eq. (101) is consistent with the classical Peach-Koehler force similarly to Eq. (78). For the case of a 2D triangular lattice or a 3D BCC crystal, the velocity takes the form with a mobility M = 2β/(φ 2 0 |b d | 2 ). The last term in Eqs. (101)-(102) accounts for the contribution from the compositionally generated stress, as a result of the compositional strain (∼ αψ) arising from local concentration variations, i.e. from solute preferential segregation (Cottrell atmospheres) around defects. The stress field may be written as with neglecting higher order terms in the last approximation obtained with η m = φ m e −iGm·u [13]. Results predicted by these equations are the deflection of dislocation glide paths, the variation of climb speed and direction, and the change or prevention of defect annihilation [13]. Simulations exploiting the FEM approach outlined in Sec. 3.3 also enable the advanced description of these effects in three dimensions, in particular for small-angle grain boundaries [13].

Multi-phase systems
Most of the APFC literature focuses on systems with a single solid phase. In a seminal work by Kubstrup et al [151], studying pinning effects between different phases, namely crystalline systems having triangular/hexagonal and square lattices, a construction has been proposed handling variable phases through a single density expansion. Extending this idea, in Ref. [58] an ansatz for the atomic density has been proposed to include more symmetries at once with {η j } and {χ m } representing different set of amplitudes associated to reciprocal lattice vectors G j and Q m , respectively. These two sets were chosen to account for the first and second modes necessary for reproducing triangular and square symmetry together, namely corresponding to J = 6 and M = 6 amplitudes. However, they can be arranged differently among the two sums, and, importantly, a reduced set of amplitudes can be exploited (see specific choices of G j and Q m in Ref. [58]). Amplitude equations would simply follow from the general equations reported in Sec. 2.3. Simulations performed with this approach, combined with the formulation illustrated in Sec. 2.4 for the excess term, showed the ability to study solidification, coarsening, peritectic growth, and the emergence of the second square phase from grain boundaries and triple junctions in a triangular polycrystalline system. See an example in Fig. 9. So far, this has been shown only for the lattice symmetry mentioned above in two dimensions. The same applies to extensions of the APFC to account for additional degrees of complexity in the crystal structure, such as for the amplitude expansion of the so-called anisotropic PFC model [124,152].

Heteroepitaxial growth
An ideal application of the APFC model is heteroepitaxial growth, where a substrate provides a single crystallographic basis for layers growing on top. In such processes, the growing film typically has similar crystal symmetry and lattice constant. The amplitudes vary on long length scales for these systems, so a relatively large computational grid spacing can be used. In this context, the large angle issue discussed in Sec. 5.1 is not present. Therefore, this would be an ideal application for using an adaptive mesh since the amplitudes in many cases vary on very large length scales. To the authors' knowledge this has not been done to date. Nevertheless, even uniform lattices can be used to study relatively large systems. An example application is a single or small number of mismatched layers grown on a substrate. The mismatch leads to interesting strain-induced Moiré patterns that have been observed in experimental systems [153][154][155]. In these cases, it is possible to model the film as a single two-dimensional layer with amplitudes. To the authors' knowledge, the largest APFC simulation of such systems was on the study of Moiré in graphene films in which the large simulation size was 19.6 µm × 34.0 µm which corresponded to roughly twenty-five billion carbon atoms. Some sample works are reviewed in the next subsection. Similarly, the amplitude expansion can also effectively be used to study the growth of many layers in two and three dimensions, i.e. to examine the Asaro-Tiller-Grinfeld [156][157][158] instability and the subsequent nucleation of dislocations. This aspect will be also illustrated in the following. This section shows the APFC model in an applied context, reproducing experimental results and outlining general properties of mismatched, multilayered systems.
6.5.1. Ultrathin films: strain induced ordering When a monolayer (or several layers) of one material are grown on a substrate, the lattice mismatch can lead to interesting strain induced patterns [159,160] and the APFC model is ideally suited to model such Figure 10. Ordering of a triangular (honeycomb) lattice on a substrate with a triangular array of potential maxima is depicted as green dots. In (a) the red, blue, pink, orange and purple dots correspond to 1 × 1 (e.g., Cu/Ru(0001) or Cu/Pd(111)), 2 × 2 (e.g., O/N(111)), (111)) and ( √ 7 × √ 7) R19.1 • respectively. In (b) the pink, red and blue atoms correspond to 1 × 1 (e.g., graphene/Cu(111)), 2 × 2 and ( patterns [28,[161][162][163][164][165]. Their nature depends on the misfit strain, ε m = (a s − a f )/a s , where a s and a f are the substrate and film lattice constants, the relative crystal symmetry of the layer/substrate system and the film/substrate coupling strength. For example, when layers of Cu are grown on a Ru(0001) substrate, the substrate potential provides a triangular (honeycomb) array of potential maxima (minima) for the Cu atoms. Since the lattice constants of Cu and Ru(0001) are similar (ε m = 5.5%), a 1×1 ordering occurs as depicted by the red dots in Fig. 10(a). For larger mismatches other orders can occur as shown in this Fig. 10 for the ordering of triangular film on a triangular substrate (TT) in (a) and a honeycomb film on a triangular substrate (HT) in (b). By symmetry a (TT) system is equivalent to a (HH) system and a (HT) system is equivalent to a (TH) system. These patterns can be characterized by two integers (k, j) or equivalently a length and angle (L, θ) as depicted in Fig. 10(a). The relationship between them is L = ja s x ( (2k + 1) 2 + 3)/2 and tan θ = √ 3/(2k + 1). In Fig. 10(a) the 1 × 1 state could occupy two equivalent separate sublattices, while in (b) this state has only one sublattice. In general, the degeneracy (N S = number of equivalent sublattices) is given by, for the TT system and half of Eq. (106) for the HT system. Figure 11(a) illustrates the different sublattices for a TT √ 3 × √ 3 R30 • system. The nature of the patterns that form depend on the degeneracy of sublattices, N , the mismatch strain, ε m , and the strength of the coupling, V 0 , between the film and substrate. In the limit V 0 = 0, a 2D Moiré pattern forms in terms of a honeycomb array of commensurate regions bounded by a triangular network of domain walls for the TT system, with length scale λ = a f /ε m . This is illustrated in Fig. 11(b) for a 1 × 1 system with a mismatch consistent with a Cu/Ru(0001). As V 0 increases, the commensurate regions increase in size, and the domain walls and junctions decrease in size but increase in energy. For the TT system, the displacement across a junction is larger than the displacement across a domain wall. Thus for the TT system at a certain V 0 it becomes energetically favorable to eliminate the junctions and form stripes. At even larger values of V 0 the film becomes commensurate with the substrate. A peculiar state in the TT arises for some values of (V 0 , ε m ) in between the stripe and honeycomb patterns in which the junction energy is lowered by twisting the domain walls and moving the junction to a lower energy location. Sample patterns for the TT system are shown in Figs. 12(a), (b) and (c). In the case of the 1 × 1 the junction energy is so high that it can create dislocation pairs and lead to zig-zag type patterns [164,165].
The HT system is considerably different since the domain wall energy is higher than the junction energy and of course the symmetry is different. At very low V 0 , a triangular network of commensurate regions forms. At a V 0 much higher than in the TT case, a stripe phase emerges. At a slightly larger V 0 , the commensurate state appears. There appears to be no equivalent twisted state in this system. Sample stripe and triangular patterns are shown if Fig. 12(e) and (f).
To model these patterns within a PFC approach and corresponding APFC it useful to consider adding an additional coupling term, F c , to the free energy functional given in Eq. (1) of the form, where V 0 is the coupling strength, the summation is over lowest order modes needed to reconstruct the symmetry of the substrate and G s m corresponds to the reciprocal lattice vectors of the substrate (which will have a different magnitude that the film). The coupling factor n j(k+1) is needed since orders greater that 1 × 1, a coupling V n would give no contribution in the amplitude expansion, since V and n would have different lattice spacings. In principle, higher order harmonics of V (or n) could be included, even though this would lead more computational expensive models. In the amplitude expansion this term leads to a coupling term F c η , of the form where D kj = ((k + 1)j)!/((kj)!j!). This term would be added to the free energy given in Eq. (30) for a triangular two-dimensional system. In addition, to account for the misfit strain, the operator G m that enters Eq. (21) becomes where α = 1 − ε m . Insight into the model can be obtained in the small deformation (u) limit, η m = φe −iGm·u . The total free energy function reduces to a two dimensional Sine-Gordon model, i.e., where C 11 = 9Aφ 2 and C 44 = C 12 = 3Aφ 2 . Unfortunately this is difficult to solve for the boundary condition of a two dimensional triangular pattern. In one dimension this reduces to a Sine-Gordon model that can be solved exactly [166]. In this model the stripe to commensurate state transition occurs when where P is a measure of the potential between the film and substrate and Ka 2 is a measure of the elastic energy in the film. These parameters are given by Figure 12. Sample patterns and phase diagrams for √ 3 × √ 3 R30 • system for TT (a)-(d) and HT (e)-(g) systems. For the TT system, the stripe, twisted honeycomb and honeycomb patterns are illustrated in (a), (b) and (c) respectively, and the phase diagram is shown in (d). Stripe and triangular patterns for the HT system are shown in (e) and (f) respectively and (g) shows the HT phase diagram. Each color in the patterns corresponds to a different sublattice. In (d) and (g) the dashed line is the analytic prediction for the stripe/commensurate transition given by Eq. (114). The figures were reconstructed from [161]. and Details of these calculations can be found in Elder et al [161]. The full phase diagram as a function of ε m and the ratio of potential/elastic energy, P/Ka 2 , can be obtained through numerical simulation. Sample phase diagrams are given for the √ 3 × √ 3 R30 • system for the TT and HT cases in Figs. 12(d) and (g) respectively. As can been seen in these figures for small ε m , the analytic predictions (this is true for all (k, j) systems) for the stripe/commensurate transition are quite accurate and very good for the HT case for all ε m .
An interesting comparison with experiments is the Cu layers on a Ru(0001) substrate which is a 1 × 1 TT system. In this case, varying the number of Cu layers increases the film's elastic energy and the potential between the substrate and film. Essentially, adding more layers corresponds to reducing the ratio P/Ka 2 . One layer forms a completely commensurate state, two layers form a striped state, three layers form a twisted honeycomb (or zig-zag state), and four layers form a honeycomb state. To compare with the non-equilibrium patterns observed in experiments, simulations starting from random fluctuations were conducted. The comparison of the experiments and simulations depicted in Fig. 13(a)-(c) shows a very good agreement for various patterns. In another experiment by Schmid et al [160] patterns of partially filled layers are reported. These patterns are remarkably similar to simulations of nonequilibrium patterns observed with the APFC model in the commensurate state as shown in Fig. 13(d).
Studies of the HT 1×1 lead to a phase diagram similar to that shown for the √ 3× √ 3 in Fig. 12. To compare with experiments, density functional theory (DFT) calculations were conducted by Smirman et al [28] to calculate the value of the dimensionless quantity P/Ka 2 for various 1×1 film/substrate systems. The phase diagram accurately predicted commensurate state for twenty-five system mostly corresponding to films consisting of monolayers of InN or GaN on various substrates. In addition, the phase diagram accurately predicted a commensurate state for graphene (G) on N, and triangular patterns for G on Cu, Pd, Pt, Al, Ag, and Au. Work was also conducted to predict the wavelength of the patterns as a function of misorientation with respect to the substrate in G/Cu (111) and G/Pt(111) systems. In the absence of coupling two dimensional patterns arise with wavelength λ = a f / ε 2 m + 2(1 − ε m )(1 − cos(θ)), where θ is the misorientation angle. The study showed that as the coupling increases, the wavelength increases and interestingly the lowest energy states were not at zero degree misorientation (0.88 • and 3.22 • for G/Cu (111) and G/Pt(111) respectively), which is unfortunately difficult to measure experimentally. However, the predicted wavelengths were consistent with the experiments of Marino et al [153] for G/Cu (111).
Other predictions of the APFC model involve the influence of defects and edges on pattern formation in the √ 3× √ 3 R30 • which corresponds to systems such as Xe/Pt (111) or Xe and Kr on graphite.
6.5.2. Epitaxial growth: island formation and defect nucleation When a material is grown epitaxially on a substrate with a mismatch strain, ε m , the film will tend to buckle and form islands or mounds as it grows due to the so-called linear Asaro-Tiller-Grinfeld (ATG) instability [156][157][158]. Recall that the APFC model is ideal for examining these phenomena, featuring relatively uniform amplitudes suited for adaptive meshing. In addition, it is possible to reduce the study of an ATG instability in a 2D film to a 1D problem [167,168]. Consider expanding about the strained film such that η m = η m e −iδqm·r where δq m is responsible for the mismatch strain imposed by the substrate. For a triangular lattice with a strain imposed in the x direction (y being the growth direction) δq 1 · r = −δ x x − δ y y/2, δq 2 · r = δ y y, δq 3 · r = δ x x − δ y y/2, δ x = √ 3/2ε m and δ y is determined by lattice relaxation. The strained amplitudes can now be expanded about a one dimensional profile, η 0 j (y) as follows and similarly for the average density about n 0 o (y) The profiles η 0 j (y) and n 0 o (y) must be determined numerically. The linearized equation of motion for the perturbed quantitiesη j andn o are quite complex but are easily solved numerically to obtain a dispersion relation (ω(q x )) for the position of the liquid/solid front, i.e., the results can be fit to the form |η j |,n o ∼ e ωt . Dispersion relations are shown in the inset of Fig. 14(a). Various analytic studies have lead to different forms of the dispersion relation depending on what physical mechanisms are included. Surface diffusion leads to ω ≈ α 3 q 3 x − α 4 q 4 x [158,169,170], wetting to ω = −α 2 q 2 x + α 3 q 3 x − α 4 q 4 x [171,172], evaporation-condensation to ω = α 1 q x − α 2 q 2 x [173,174] and bulk diffusion to ω = α 2 q 2 x − α 3 q 3 x [175]. In the APFC simulations, ω can be fit to a fourth order polynomial in q x however none of the fits are consistent with any of the prior results. This is due to the fact that the APFC model cannot separate each of the mechanisms individually.
From these studies the most unstable q x , Q * , can be extracted as a function of misfit strain and interface width (W ) as shown in Fig. 14(a). The width, in the notation of Eq. (4), was altered through the variable B x since W ∼ B x /|∆B 0 | [45]. For small values of ε m it was found that Q * ∼ ε 2 m and for larger values Q * ∼ ε m for all interface widths. ATG theory gives Q * ≈ (E/γ)ε 2 m , where E = B x φ 2 /2 is Young's modulus, φ Figure 14. (a) most unstable wavevector (Q * ) is shown as a function of misfit strain (ε m ) for various interface widths. In the inset dispersion relations are shown for ε m = 4 % (red) and 3 % (blue). (b) the Q * and ε m are rescaled to give rise to a universal curve as described in the text. In the insetQ is shown as a function ofε 2 m . Details of the calculations can be found in reference [168]. Reconstructed from [167,168].
is the magnitude of the amplitudes in equilibrium, γ is the surface energy which can be calculated numerically. The numerical results fit the small ε m to Q * = 4Eε 2 m /3γ. The linear behavior at large ε m can be understood by considering the wavelength at which the insertion of a dislocation would lead to perfect relaxation (i.e., the addition or subtraction of a lattice point every λ returns the lattice constant of the film to its equilibrium value). This occurs when Q * = 2π/λ = q x |ε m |. It is interesting to note that this linear relationship was observed in experiments on SiGe/Si(001) growth [176,177] although other explanations may exist as this is a binary system [178].
The continuum (ATG) calculation fails when the most unstable wavelength (2π/Q * ) becomes comparable with the interfacial thickness. If one supposes that the crossover occurs at ε c m when 4Eε c m /3γ = q x ε c m then ε c m = 3γq x /4E and Q c = 3γq 2 x /4E. Defining the scaled quantitiesε m = ε m /ε c mQ = Q * /Q c gives rise to the universal behavior shown in Fig. 14 (b). That is, the relationship betweenε m andQ is independent of the interfacial thickness. It was found numerically that 1/Q c ∼ B x ∼ W 2 .
An APFC study of the growth of islands of one material on a ribbon of another was conducted by Elder et al [6,150]. Several experiments [179][180][181] had to be undertaken to examine whether the growth of islands (or quantum dots) on thin ribbons may be exploited for better control of island sizes and correlations. When an island of one material grows on an island of another material, the misfit strain will eventually lead to the nucleation of dislocation at the island/film/vapour junction. On very thin ribbons, the strain in the island can be somewhat reduced by bending the ribbons, leading the possibility of growing larger defect-free islands. An example is shown in Fig. 15. Figures  Figure 15. In this figure the magnitude of the sum of the amplitudes is shown for an island of one material grown on another. In (a)-(c) the time evolution of one island is shown. Similarly in (d)-(f) an island growth is illustrated for a thicker ribbon. In (g)-(l) the time evolution of island growth and nucleation is shown. In (a)-(f) a flux of material only came from the top, while in (g)-(l) it came from both sides of the ribbon. Reconstructed from [6].
(a)-(c) and (d)-(f) show the growth of an island for two different ribbon thicknesses. In (c) and (f), the final island size (L f ) at which dislocations appear indicates that L f is larger for the thinner ribbons. Depending on conditions it was shown in reference [150] that decreasing the ribbon size could almost double L f . Another interesting feature emerges when the island starts to grow. It bends the ribbon such that preferential regions for island nucleation appear on the other side near the triple junctions, leading to correlated growth as shown in Fig. 15(g)-(l). This correlation could potentially be exploited to create uniform arrays of islands.
In summary, the binary and pure APFC models provide an excellent platform for studying heteroepitaxial growth. Coupled with adaptive mesh schemes as illustrated in Sec. 3, very large simulations should be possible in both two and three dimensions.

Conclusions and outlook
In recent years, bridging-scale modeling has become crucial to comprehensively investigate crystalline systems, explore macroscopic effects of microscopic details, and unveil general properties and behaviors for further scale-specific characterizations. Here, an overview is provided of the model(s) obtained through the amplitude expansion of the phase-field crystal (APFC), which combines the description of crystals on relatively large (diffusive) time scales, conveyed by the PFC model [1,7,33], with a spatial coarsegraining. The concepts underlying its derivation have been illustrated, focusing on practical aspects such as explicit formulas, generalizations, and examples, along with presenting different formulations.
Computational aspects have also been outlined. The fields (amplitudes) to solve for within the APFC model are suited for inhomogeneous spatial discretizations, a feature that motivated its development in the first place [23]. Recently, a few optimized methods have been developed to allow for large-scale calculations and, in particular, paving the way for extensive three-dimensional calculations.
The APFC model emerges as one of a kind among mesoscale approaches: it handles the description of crystalline systems through slowly varying continuous fields, so without resolving atoms, but retains details of the crystal structure such as anisotropies and lattice defects. Namely, it merges different aspects addressed by micro-and macroscopic approaches within a single model rather than coupling models working at different time-and length scales (like other remarkable approaches as, e.g., the quasi-continuum approach [182,183]). Among its key aspects, special attention has been given to the mesoscale description of elasticity and plasticity, being the primary goal of many coarse-grained descriptions (as the phase-field crystal itself [1,7]). As a pivotal example, the elastic field generated by dislocations within the APFC model matches classical continuous descriptions and encodes a core regularization related to the lattice parameter. Moreover, it is expected to be affected by lattice symmetry and encodes nonlinearities. Amplitudes also allow for characterizing plasticity and defect dynamics. This description can be exploited within the broader context of PFC models as amplitudes fully characterize deformations therein [75]. Like every other model, APFC has its range of applicability, strengths, and weaknesses. One weakness is the ability to accurately predict the precise structure of atomic-scale structures such as dislocations and interfaces, similar to the drawbacks of traditional phase field models. However, it may be employed to investigate long-range effects for such systems, and extensions have been provided to improve the mesoscale descriptions with respect to the standard formulation (see, e.g., the control of energies for defects and interfaces and the modeling of Peierls barriers). Like PFC, the variational, overdamped formulation of the APFC model conveys a lack of separation among different timescales, affecting the competition among diffusion mechanisms and elastic relaxation. This issue, however, has been solved by a few different extensions, which are expected to become the standard approaches for phenomena when the separation of timescales is relevant. The most critical aspect for applications of the APFC model remains the limitation to small rotations with respect to a reference crystal orientation (see the problem of beats [66,73,74]). It prevents the thorough investigation of high-angle grain boundaries and polycrystalline systems. Therefore, providing a solution for this issue is a crucial challenge for achieving a general mesoscale description of crystals. To date, this aspect has been only partially addressed through a covariant formulation with respect to rotation of the crystals, which still needs to be assessed for the description of elasticity and plasticity and its compatibility with other extensions.
It is worth mentioning that in light of the limitation(s) mentioned above, the currently available APFC models should be considered valid for relatively small deformation and rotation only, de-facto for every crystalline system where defects as dislocations can be described as separated objects. However, systems featuring such conditions are common, widely studied, and exploited in several technology-relevant applications, such as single crystals, alloys, and homo-/ heteroepitaxial systems, besides small angle-grain boundaries. The overview and discussion of the main applications addressed so far in the literature illustrate this aspect.
In conclusion, this review has attempted to collect the basics and the recent developments of the APFC model. While it has been used to study several physical phenomena, its potential still has not been fully exploited. Potential applications include the investigation of three-dimensional mesoscale tracking of defects and interfaces (e.g., for heteroepitaxial systems). Moreover, besides the challenges already mentioned above, a few aspects can be identified which will improve the approach further: i) direct connections with advanced continuum theory for elasticity and plasticity, closing the gap with methods such as dislocation dynamics; ii) description of complex crystal symmetries beyond simple ones to broaden the application to technologyrelevant systems; iii) extending the parametrization to include physical parameters extracted from experiments and/or other methods; iv) connections and coupling to both microscopic, fully atomistic (e.g., PFC or Molecular Dynamics) and macroscopic (e.g., phase-field, continuum elasticity) models; v) extended boundary conditions to enable investigations beyond bulk-like systems and simple geometries; vi) further development of numerical methods, keeping up with state-of-the-art numerical techniques.