Explicit temperature coupling in phase-field crystal models of solidification

We present a phase-field crystal (PFC) model for solidification that accounts for thermal transport and a temperature-dependent lattice parameter. Elasticity effects are characterized through the continuous elastic field computed from the microscopic density field. We showcase the model capabilities via selected numerical investigations which focus on the prototypical growth of two-dimensional crystals from the melt, resulting in faceted shapes and dendrites. This work sets the grounds for a comprehensive mesoscale model of solidification including thermal expansion.


Introduction
Solidification of crystalline materials is an ubiquitous phenomenon in nature and technology. It yields shapes and patterns which depend on out-of-equilibrium growth conditions and the interplay of different instabilities [1]. At the same time, the control of the solidification process is key in several technological applications, from conventional ones [2,3] up to bottom-up approaches exploiting self-assembly [4,5]. The outcome of crystal growth results from several different contributions: crystal seeds grow after nucleation or on pre-existing crystal phases/seeds, while later stages are affected by capillarity, elasticity, plasticity, as well as various kinetic effects [6][7][8][9][10]. When considering solidification, heat transfer in the region surrounding the solid phase also plays an important role [6,[8][9][10]. From a microscopic point of view, the arrangement of atoms in a periodic lattice typically emerges in anisotropic behaviors such as faceting, as well as in the nucleation and motion of defects, which are tightly related to crystallographic directions.
The modeling of solidification has been addressed by different approaches. Microscopic approaches are suitable for evaluating lattice-dependent features such as anisotropies and defect structures [7]. However, the growth of crystals involves long time scales, typically not accessible within these methods. Macroscopic approaches, which cope with large systems and long timescales, proved successful in describing the main features of crystal growth via both front tracking and phase-field methods [11][12][13][14][15][16][17]. However, they usually lack a direct connection to the lattice symmetry and microscopic features in general. Lattice-dependent effects can then be described only partially and are included mainly through parameters and additional functions, e.g. as anisotropic interface energies [18,19].
The so-called phase-field crystal (PFC) model [20][21][22][23] emerged as a prominent approach to describe crystal structures at large (diffusive) timescales through a continuous, periodic order parameter representing the atomic density. It describes solidification and crystal growth, including capillarity, elasticity as well as nucleation and motion of defects. As a result, it reproduces the main phenomenology of crystal growth in two-and three-dimensions [23][24][25]. Also, it allows for a self-consistent description of anisotropies resulting from the lattice structure [26,27]. Therefore, it represents a suitable framework for the development of a comprehensive model of crystal growth. Only recently, however, a temperature field and the related heat flux have been explicitly considered in PFC models [28,29]. Similarly, while its first formulation already includes elasticity effects, recent studies shed light on its connections to continuum mechanics [30][31][32][33][34][35][36], allowing for extended descriptions of the crystal lattice and its deformation.
This paper presents a PFC model of crystal growth accounting for temperature effects through the coupling with thermal transport [29], and including a temperaturedependent lattice parameter by a suitable modification of the PFC free energy [36]. The resulting model is illustrated in section 2. The equations are solved numerically by a simple but efficient approach based on the Fourier pseudo-spectral method as shown in section 3. Section 4 reports the description of the continuous elastic field exploited to assess the proposed model. In section 5 the model is then showcased through selected numerical simulations. They reproduce the prototypical growth of twodimensional crystals from a melt, resulting in faceted shapes and dendrites. The effects of temperature described by the model are illustrated. Finally, we draw conclusions and perspectives in section 6.

Model
The PFC model is based on a Swift-Hohenberg free energy functional [20,21,23], which can be written as with ψ ≡ ψ(r, t) a scalar order parameter related to the atomic number density, Ω ∈ R n its domain of definition (n = 2 in this work), q 0 a length scale enforcing a periodicity of a 0 ≈ 1/q 0 and a set of material parameters λ, κ > 0. Coupled with appropriate boundary-and initial-conditions the time evolution of ψ is described via a conservative (H −1 ) gradient flow of F , namely It is worth mentioning that different forms of F are used in literature as, e.g., in the originally proposed formulation [20] including q 0 in the differential operator ψ(q 2 0 +∇ 2 ) 2 ψ or (q 2 0 ψ +∇ 2 ψ) 2 only, a parameter = κ−λ playing the role of the quenching depth, and eventually additional coefficients for other terms in ψ, such as τ ψ 3 + νψ 4 (with τ = 0 and ν = 1 in [20]). Such formulations, however, are equivalent and lead to the same expression for δF/δψ. In [29] an extension of this PFC model including heat-transfer through a temperature field was proposed. Therein, the (dimensionless) Helmholtz free energy functional with q 0 = 1, here normalized with respect to the temperature, reads with T ≡ T (r, t) the dimensionless temperature and T = 1 at the melting point. ϑ, γ > 0 are additional parameters. The dynamics is then given by the coupled evolution of T and ψ, reading with M > 0 a parameter corresponding to the thermal diffusivity, assumed to be constant. This model is here extended by including a variable lattice parameter and connecting it to the temperature field, aiming at describing thermal expansion in the crystal phase. Let us consider a crystalline solid with lattice parameter a 0 . Expansion due to temperature can be described through the thermal expansion coefficient α th = L −1 (dL/dT ) withL a characteristic (reference) length of measurement and L a linear dimension in the solid. Assuming weak dependence of α th on T , a T = a 0 [1+α th (T −T 0 )] with a T and a 0 the lattice spacing at T and T 0 , respectively. This equation defines a relative deformation, i.e. a thermal strain ε th (T, T 0 ) = (a T − a 0 )/a 0 = α th (T − T 0 ). Following the arguments in [36], this information may be encoded in the PFC free energy by noticing that (a T − a 0 )/a 0 = (q 0 /q T ) − 1 and thus q T = β(T, T 0 )q 0 with β ≡ β(T, T 0 ) = 1/(1 + ε th (T, T 0 )). Therefore, by setting q 0 = 1 we may derive from equations (1) and (3) the following free energy functional: This leads to the following evolution equations defining our model: Note that the PFC heat-flux model, equations (3)-(4), without temperature-dependent lattice parameter [29] is recovered for β = 1. On the other hand, by neglecting the heat flux with a constant and unitary temperature field, and setting a constant ε th , the evolution equation for the density ψ in the presence of an eigenstrain ε * = ε th as introduced in reference [36] is obtained. Therefore, our resulting model, equation (6), realizes a general combination of these two models.
In this work we consider honeycomb and triangular two-dimensional crystals. In a one-mode approximation, they are well described by with i the imaginary unit, ψ ≡ ψ(r, t) the local average density, η n ≡ η n (r, t) amplitudes function which are real for relaxed crystals, η n = φ, and reciprocal lattice vectors According to energy minimization, the real amplitude in the relaxed solid phase results with A the global average density. For A > 0.5 (A < 0.5) the amplitudes can be identified with φ − (φ + ) resulting in honeycomb and triangular structure, respectively (with similar energetics). Equation (8) is exploited to initialize ψ in the solid phase, whereas the density of the pure liquid phase is initialized as ψ = A and φ = 0. Also the temperature field is initialized as spatially homogeneous T = T 0 . Parameters used in the simulations are reported in table 1.

Numerical methods
Numerical solutions of the partial differential equations (6) are computed using a firstorder linear IMEX scheme for the time discretization in combination with a Fourier pseudo-spectral method for the spatial discretization, enforcing periodic boundary conditions. The numerical simulations are performed on a uniform grid with an element size of dx = dy = 0.78125, resulting in a resolution of ≈ 75 grid points per unit cell. The Fourier transforms are performed with an algorithm based on the FFTW2 library. A constant time step size ∆t for the IMEX scheme is considered, and chosen to ensure numerical convergence. Additionally, we exploit the numerical time-stabilization routine presented in [37]. This approach features a convex-concave splitting of the free energy controlled by a parameter C, combined with an IMEX scheme (C-IMEX). For every set of model parameters, we determine optimal C by minimizing the difference of the error (least square of the difference) in the free energy decay obtained with the C-IMEX scheme with respect to a numerical reference solution (IMEX scheme with small timestep). Optimal errors are found by considering a minor rescaling of the time scale by a factor K ≈ 1, which is found to compensate for quantitative effects on the time scale observed for the C-IMEX scheme [37]. A convergence study of the considered numerical methods is reported in figure 1, for a specific set of model parameters leading to the growth of a crystal in domain Ω = [−28p, 28p] 2 and for time steps ∆t ∈ [10 −2 , 1]. This figure shows the residua of ψ and T , R ψ and R T , respectively, evaluated as the integral over the domain of the squared difference from a numerical reference solution for different ∆t. The IMEX (C = 0) and C-IMEX (C = −0.149) schemes are compared, with a reference solution corresponding to the C-IMEX scheme with ∆t = 0.01. As expected, the schemes converge linearly for a decreasing time step size ∆t. However, for a fixed ∆t, the C-IMEX approach gives two orders of magnitude smaller error than the IMEX approach for ψ and a two and a half orders of magnitude smaller error for T . By fixing an accuracy instead, the C-IMEX approach allows for two-order of magnitude larger timesteps compared to the IMEX approach. Convergence studies similar to figure 1 have been performed for all  Table 1. Numerical parameters (C, K, ∆t) for all the simulations reported in this paper. We set M/γ = ϑ/γ = 1 except for γ = 0 (decoupled equations in (6)), for which we set M = ϑ = 0. Empty table entries read as the row above.
the simulations reported in this work. Table 1 summarizes the numerical parameters (C, K, ∆t), together with other model parameters for all simulations in this paper.

Evaluation of elastic fields
In this section we address the evaluation of the elastic field to assess the PFC formulation including a variable lattice spacing as introduced in section 2. Elasticity in PFC models can be characterized by small perturbations of the equilibrium density ψ. Therefore, it features variations of ψ over a length scale significantly larger than the lattice spacing, and a coarse-graining of ψ is usually considered in this context [35,38]. In particular, one may consider an expansion of the periodic density as in equation (7) with complex amplitudes encoding lattice distortion over a displacement field u. Indeed, setting (7) with η n (r, t) = φ(r, t)e −iq n ·u(r) . This ansatz can be also exploited to derive a coarse-grained PFC model, namely the amplitude expansion of the PFC (APFC), where η n are the variables to solve for [39][40][41]. Within the PFC model, η n as well as ψ can be computed from a demodulation of ψ [38,42]: with F the Fourier transform, F (−1) the inverse Fourier transform, and implying Einstein summation convention. The demodulation (9) includes a coarse-graining procedure over one unit cell UC = [0, Following recent works [33,36,38,43], one may obtain the mechanical stress tensor σ ij ≡ σ ij (r, t) without isotropic pressure terms from the density field: with · a spatial average over the unit cell, which may be explicitly performed in reciprocal space through the smoothing kernel similarly to equation (9) for ψ. Equation (10) generally defines σ ij up to a phase-dependent divergence-free contribution, which can be removed by considering ψ as expressed through an onemode approximation as in equation (7) [38]. A coarse-grained description of the crystal is also achieved through amplitudes from equation (9). Assuming η n to be constant over the length of UC (similarly to the basic assumption of the APFC model [41]) and denoting the complex conjugate of η n as η * n , we obtain Equation (10) with an one mode approximation for ψ and equation (11) are found to deliver almost identical elastic fields. As will be shown below, this approach allows to characterize elastic properties also in systems featuring solid-liquid coexistence.
In reference [36] a benchmark problem has been reported for the APFC model. It consists of the simulation of an elastic inclusion through an eigenstrain formulation that enforces a spatially-dependent lattice parameter through a quantity similar to β in equation (5). It corresponds to the so-called Eshelby problem [44][45][46]. Here we consider the same setting to assess this elastic inclusion problem within the microscopic PFC model (rather than the coarse-grained APFC model). We set a spatially dependent, but constant in time, lattice parameter √ 3p, 25 √ 3p] for a triangular lattice and we let the system relax until a steady state is reached. Figure 2 provides a comparison of the mechanical stress components computed from the PFC simulation through equation (11) and the analytic solution for an inclusion in an infinite medium [46,47]. For comparison, the result obtained by an APFC simulation [36] is also reported. A good agreement between the analytic solution and the simulations is obtained. A few differences may be recognized and ascribed to different effects similarly to reference [36]. In brief, the considered setting models the inclusion problem inside a domain of finite size and with periodic boundary conditions, while the Eshelby solution is derived for infinite systems. Also, elastic properties of PFC models account for non-linear effects, which are absent in the analytic solution. Furthermore, the approximation of the inclusion's boundary as a diffuse interface achieved through equation (12) inherently produces a smooth field that differs from the sharp transition implied in the analytic solution. This smoothing and the decay of the elastic field away from the inclusion match the results obtained with the APFC model with the same parameters. A slightly different stress field in the inclusion is obtained, that can be ascribed to a fixed average density assumption in the considered APFC simulation.

Model features and parameter study
In this section, we illustrate the capabilities of the proposed model through numerical experiments exploiting the method described in section 3. In particular, we provide an overview of its features focusing on the dendritic solidification regime. Figure 3(a)-(b) illustrates the numerical solution of equation (6) at a representative stage (t = 3.5 · 10 4 ) for A = 0.849 (A > 0.5, honeycomb), with and without including a temperaturedependent lattice spacing (namely with α th = 0 and α th = 1.0). Focusing first on the morphologies of the crystalline domain, we note that the initial crystal seed grows in both cases and develops an anisotropic shape. It corresponds to a six-fold dendrite which reflects the underlying triangular symmetry of the crystal lattice. This evidence extends the preliminary numerical results reported in reference [29], showing that the model including heat flux allows for dendritic growth with an extended set of parameters. Moreover, differences are observed among the simulations with α th = 0 and α th = 1.0. The latter case leads to growing arms of the dendrite, which are more rounded than the former. As shown in figure 3(c) reporting a comparison of the length of dendritic arms over time D ≡ D (t) in both cases, the velocity of the tip is also smaller for α th = 1.0 (namely D /t ≈ 1.85·10 −2 and D /t ≈ 1.50·10 −2 for α th = 0 and α th = 1.0, respectively). However, in both cases, after a first initial transient, a constant tip velocity is observed, consistently with a classical result for dendritic growth [10]. Figure 3(c) also includes the time evolution of F β , showing that the energy decreases for both the considered settings.
For the chosen parameters, the growing crystal has a local average density lower than the liquid phase and, as dictated by the equation for ∂ t T , a lower temperature with respect to T 0 (see figure 3(a)-(b)). This behavior is observed, e.g., for materials such as hexagonal ice which is well represented by a honeycomb structure, but also diamond-cubic silicon, rhombohedral bismuth, and orthorhombic gallium [48]. At the liquid-solid interface, the temperature increases above T 0 , wheres away from the growing seed where ψ = A, we have T ≈ T 0 . These variations can be interpreted in terms of a Gibbs-Thompson effect [49,50], commonly observed at phase boundaries, combined with the conservation laws. A region with a slight material depletion surrounding the T > T 0 region is observed for α = 1.0, which can be ascribed to a material transfer towards the interface triggered by the reduced lattice spacing in the solid phase. In turn, larger temperature gradients are present, which enforce the more isotropic shapes. The more typical behavior of larger (local) average density and T > T 0 in the crystal phase can be straightforwardly obtained by setting a global average density A → 1 − A (i.e., for the corresponding triangular phase, as can be readily shown from the equations defining the model). These cases lead to identical energetics and, in turn, morphologies and fields for the α th = 0 case. However, for α th = 1.0, either a lattice expansion or compression is enforced for settings with symmetric global average densities with respect to A = 0.5. As a result, differences are observed: for A < 0.5 no slight material depletion is observed from the liquid phase with instead a larger region with an increased average density, temperature gradients closely resembling the α th = 0 case and, in turn, similar morphologies (not shown), in agreement with the arguments reported above. This symmetry breaking will be exploited in the following section to characterize the thermal stress. More insights into the morphological evolution are obtained through a parameter study for the same simulation setup considered above. In particular, to shed light on temperature-dependent effects, we vary T 0 , γ, α th separately in the following range: We analyze D as well as the temperature range, ∆T (t) ≡ max x∈Ω T (x, t) − min x∈Ω T (x, t), and the morphologies at a given time (t = 3.5 · 10 4 ), as illustrated in figure 4. Deacreasing T 0 or increasing γ leads to a decrease in ∆T and D at a given time, and vice versa. Also, small changes in the shapes are obtained as illustrated in insets of figure 4(a)-(b) comparing representative morphologies. We conclude that the initial temperature, T 0 , and the coupling strength between equations (6), γ, mainly affect the growth rate. In contrast, a large thermal expansion coefficient, α th , leads to more significant differences in terms of morphologies consistently with the result in figure 3, as well as a decrease of D and increase of ∆T . Still, the growth dynamic results mostly unaffected for α th < 0.1. While kept fixed in the simulations discussed above, parameters entering the original PFC model, like A, λ, κ, as well as the size of the simulation domain, still determine the features of the solidification process and if the growth leads to a dendritic shape. For the sake of completeness, we illustrate in figure 5 the main expected changes for the solidification process. Varying the average density affects the anisotropy of the growing crystal/dendrite. Also, effectively changing the quenching depth through λ and κ may lead to a different growth regime, up to reaching equilibrium with a morphology corresponding to the anisotropic equilibrium crystal shape, see figure 5(a). It is worth mentioning that the growth of dendrites discussed so far may be considered in a larger domain allowing for the development of long arms and the onset of the growth of secondary arms. This is illustrated in figure 5(a)-(c).   (6) with parameter ranges as in (13). The effects of (a) T 0 , (b) γ, and (c) α th are illustrated. Insets show representative comparisons of morphologies with a reference (corresponding to the empty triangle in every panel).

Crystal growth and lattice deformation
After characterizing the model, we study the growth in a prototypical case representative of more general settings. We consider the solidification occurring in a system with two crystal seeds having different (random) orientations. We address both the cases with (α th = 1.0) and without (α th = 0) a temperature-dependent lattice spacing. Figure 6 shows the resulting density-and temperature fields at four representative stages. The dendrites exhibit shapes analogous to figure 3 owing to the same parameters.  Moreover, when the length of the arms approaches half the distance among the seeds, the temperature fronts where T > T 0 meet, and the two growing crystals begin to interact. At this stage, similarly to results reported in the literature [17,51], the growth rate drops. Interestingly, while for longer times, the dendritic arms are still separated in case α th = 0, merging is observed for the α th = 1.0 case, with the formation of grain boundaries and dislocations. In general, the merging of growing crystals depends on the phase of the periodic lattice and its deformation at the corresponding growth fronts. If they would form a relaxed crystal, i.e. if the crystal structures meeting at the growth fronts are commensurate, the process is expected to occur. Otherwise, additional energy barriers exist, which may lead to the hindering of the merging process. The evidence from figure 6 suggests that the effects induced by the temperature dependence of the lattice spacing in the specific setup considered for this simulation favor the commensurability of the crystals at the growth fronts.
To understand the behavior observed in figure 6 and further characterize the effects described by the model, we focus on an idealized case, namely the growth of a small seed and the interaction with its periodic images in a domain hosting an integer number of the unit cell at T 0 . The growing crystal is expected to interact with its periodic images. As shown in figure 7(a), with α th = 0, no merging occurs, and tensile stress is observed (see figure 7(c)). The latter builds during the very first stages of the growth and decreases up to a given value when the temperature equilibrates to T 0 . Note that no significant differences are observed in this regard with (γ = 0) and without (γ = 0) the coupling with the heat flux (see figure 7(c)), so this effect is present in the classical PFC model. With α th = 1.0, we consider two cases corresponding to two different initialization of the average density. We consider A 1 = A, for which T < T 0 in the crystal, and A 2 = 1 − A, for which T > T 0 in the crystal, such to have opposite T − T 0 contributions (see also discussion of figure 3). These two cases show different behaviors. The former case features lower tensile stress due to a negative contribution of the temperature-induced lattice expansion, partially compensating for the tensile stress observed during growth. As a result, merging occurs, and smaller residual stress is observed at later stages. The latter, instead, features an additional tensile-stress contribution due to the increase in the temperature, and no merging is observed. In this case, similar tensile stress is observed owing to the equilibration of the temperature at T 0 for an isolated crystal with a shape very similar to the α th = 0 case. Exploiting the simulations initialized with both A 1 and A 2 we may provide an estimation for the thermal stress. The mechanical stress for α th = 0, namely σ 0 ij may include different contributions in general, e.g., effects of an applied load and specific boundary conditions. Under similar conditions but with α th = 1.0 an additional contribution resulting from the dependence of the equilibrium lattice spacing on the temperature field, σ α ij , is expected such that σ ij = σ 0 ij + σ α ij . We point out that σ 0 ij is not necessarily equivalent for α th = 1.0 and α th = 0 as the growth dynamics is different. However, it is expected to be similar for early stages of the growth where the morphology do not differ significantly, as we consider here. An estimation of σ α ij may be then given by defining σ ij (A k ) = σ 0 ij (A k ) + σ α ij (A k ), with k = 1, 2, and σ α ij (A 2 ) ≈ −σ α ij (A 1 ) =: σ α ij , reflecting the behavior of the temperature in the solid. We may then compute ]. In figure 8 we show the mechanical stress components for the simulations initialized with average densities A 1 and A 2 , and address their splitting into σ 0 ij and σ α ij for an early stage of the simulations reported in figure 7 with α th = 1.0. Independently on the value of α th , the stress vanishes in the liquid phase, while its values changes for A 1 and A 2 , owing to the different temperatures, see Fig 8(a)-(b). More insights are given in Fig 8(c), reporting σ xx (A 1 ) and σ xx (A 2 ) along the x-axis crossing the center of the crystal as well as the computed σ 0 ij and σ α ij . σ 0 ij closely resembles the value obtained for α th = 0 case. Importantly, opposite effects are obtained for T < T 0 and T > T 0 in the solid obtained with A 1 and A 2 , respectively. A positive (negative) σ α ij is obtained for T > T 0 (T < T 0 ) in the crystalline phase, consistently with the expected thermal stress. A small region with opposite behavior at the solid-liquid interface may be ascribed to the features of the solidification process [52,53] and includes the effects induced by changes in the temperature field. It is worth noting that considering a proper separation of the timescales for elastic relaxation and diffusive dynamics, eventually controllable through parameters, would account for a complete model allowing for the exploration of different relaxation regimes [32,33,38,54]. However, the purely diffusive dynamics considered here, delivering a slow elastic relaxation as encoded in the standard PFC, is instrumental in interpreting the different quantities entering the model, which may be extended without major adaptations of the aspects discussed above.

Conclusion
We introduced a PFC model of crystal growth accounting for heat fluxes and thermal expansion. As a starting point, we considered a thermodynamically consistent approach coupling the evolution of the density field with a diffusion-reaction equation for the temperature field. Then, a parameter controlling the periodicity of the density, dependent on a temperature field has been set consistently with the classical law of thermal expansions in solids. We assessed this description by comparison with a known analytic solution of a prototypical system, namely addressing the Eshelby inclusion problem in the PFC model. Moreover, we showed the model's capabilities by numerical simulations performed with a simple but efficient numerical scheme. An overview of how the parameters entering the model control the growing morphologies is provided, focusing on the dendritic growth regime. Importantly, we characterize thermal stress effects in the system and the change of spacing induced by temperature changes. This is found to affect the merging process of growing crystals.
This work sets the ground for a comprehensive approach to crystal growth at diffusive time scales, retaining microscopic features. Perspective extensions include the modeling and simulations of the growth of faceted crystals and dendrites in open systems beyond closed ones enforced here by periodic boundary conditions and accounting for the latent heat of solidification. Like classical PFC and phase-field models, the approach presented here could be readily applied to the study of three-dimensional systems with the aid of efficient numerical methods [25,55]. However, for this purpose, the so-called amplitude expansion of the PFC model [39,41] may also represent an ideal framework as it features an additional spatial coarse-graining and may enable large-scale simulations. Additionally, considering an explicit modeling of elastic relaxation may allow for a better description of competitive time scales, in particular concerning elastic relaxation and diffusive dynamics [32,33,38,54], also in the presence of heat flux.