Phase Field Modeling of Dendritic Growth on Spherical Surfaces

The spontaneous formation of a crystal phase is one of the most common and beautiful pattern formation mechanisms in nature. Different instabilities in the crystal interface may lead to the growth of ramified structures, known as dendritic crystal growth. In this work, we use a Phase Field Model and numerical simulations to study 2D dendritic growth on curved surfaces. We show how the degree of ramification of a growing nucleus is modified by the underlying curvature of the substrate.


INTRODUCTION
The formation of a crystalline structure from an initial liquid phase (crystallization), is one of the richest pattern formation processes, with deep consequences in condensed matter physics and chemistry, material science, and even biology (Kashchiev, 2000;Kelton and Greer, 2010). For years this process has been scrutinized by means of experiments, simulations, and mathematical models, in both 3D and 2D systems.
In general, the fundamental mechanism of crystallization is given by nucleation and growth, where an initial fluctuation in the liquid spontaneously forms a small seed of the equilibrium crystal. In the simplest picture, there are only two energies associated with this seed (Kashchiev, 2000;Kelton and Greer, 2010). There is an energy gain due to the formation of a piece of the equilibrium (less energetic) structure. But there is also an energy penalty due to the formation of an interface between the crystal and the liquid. This competition leads to an activated process where only those seeds overcoming a critical size are able to grow, while the rest collapse and disappear by surface tension.
Conventional nucleation and growth commonly leads to the formation of growing compact nuclei, with typical polygonal shapes depending on the anisotropies of the surface tension. However, instabilities in the liquid-crystal interface, originated for example by the diffusion of heat needed in the formation of the crystal phase, may lead to the ramifications of the nuclei, in a process known as dendritic crystal growth (Langer, 1980).
On the other side, in the last years, condensed matter scientists have also begun to study the properties of 2D ordered phases deposited on curved substrates (Nelson, 2002;Vitelli and Nelson, 2004;Vitelli et al., 2006;Irvine et al., 2010;Tarjus et al., 2012;García et al., 2013García et al., , 2015Gómez et al., 2018). In such cases, it has been shown that the underlying curvature of the substrates modifies the structure of the phases, such that topological defects, like dislocations and disclinations, equilibrate and locate on particular regions of the substrate. This is because, on curved geometries, topological defects may help to reduce the strain energy in the ordered phases (crystal or liquid crystal phases).
In the same line, it was recently shown that the curvature of substrates can also affect the crystallization dynamics and mechanism of crystal growth on curved surfaces (Meng et al., 2014;Horsley et al., 2018). This is because curvature modifies the imbalance of volume and surface free energies of nuclei leading to crystallization. For example, in the case of nucleation on spherical substrates, it was shown that the size of a critical nucleus R c on a sphere of radius a decreases, as compared with the critical size (at same conditions) on the plane R 0 c , following the expression : Thus, it is easier to nucleate on a sphere than on a plane. This is just because on a sphere a nucleus has a bigger ratio between volume and surface energy, favoring the nucleation of the equilibrium phase. Such expression, obtained from a coarse grained phase field model, has been confirmed by molecular dynamic simulations (Law et al., 2018), and complementary theoretical calculations (Horsley et al., 2018). On curved surfaces, it was also found that curvature modifies the shape of critical and growing nuclei (Meng et al., 2014;Köhler et al., 2015Köhler et al., , 2016Paquay et al., 2017;Ma et al., 2019,?;Wang et al., 2019). This is mainly observed when the growing nuclei are unable (due to the particles' interactions) to relax the strain energy by nucleating topological defects (dislocations and disclinations). In such cases, curvature induces and elastic instability which produces the growth of ramified un-defected nuclei.
There is also an interest in understanding the features of irregular crystal growth on curved surfaces. In this sense, the dendrite growth of metal patches on colloidal particles was recently experimentally studied (Bihr et al., 2017), and it was also experimentally showed that the mechanism of freezing of soap bubbles (crystallization on spheres) may involve convection flows transporting seeds and small crystallites across the system (Ahmadi et al., 2019). In addition to that, crystal growth on the non-plane surface has proven to be relevant in understanding the structures of virus (Ganser et al., 1999;Dharmavaram et al., 2017) Here, we use numerical simulations to study the features of adiabatic crystal growth on curved surfaces. For the sake of concreteness, and due to the experimental relevance, here we will focus on crystal growth on spherical substrates. By using simple models, we will show how the curvature of spheres modifies the growth of compact and dendritic (ramified) nuclei.

MODEL
Here we study crystallization in curved geometries within the framework of the Phase Field Model approach (Provatas and Elder, 2011), by using a scalar and real order parameter φ(r, t). Points on the surface of the substrate are specified by curvilinear coordinates r = (x 1 , x 2 ). In this coordinates, the metric of the surface takes the form ds 2 = g αβ dx α dx β , where g αβ is the metric tensor (O'Neill, 1997). The order parameter φ is a measure of the degree of order in the system during the phase transition, its values ranged from 0 when the system is in the disorder (liquid) phase, to 1 when the ordered (crystal) phase has been formed.
In this model, the phase equilibria is described in terms of the total free energy F of the system expanded as a functional of φ(r, t). In general, the free energy for a mixed state system is express as (Yu et al., 2017): The term e(r) represents the magnitude of penalization in the free energy by the interphase, and f (φ, T) is the local free energy density of a homogeneous system having an order parameter φ at temperature T, and g is the determinant of the metric tensor. The dimensionless temperature T is such that T = 0 is the subcooling temperature, and T e = 1 is the two-phase equilibrium temperature. The function f has the typical shape of a double well potential with two local minima, one of them φ = 0 corresponding to the liquid, and the other in φ = 1 to the crystal phase. The competition between the desire of the system to remain in one of the bulk phase minima of f (liquid or solid) and the cost of high gradients results in a finite interface width (Warren et al., 2003). The local free energy f is represented as: where the parameter m(T) controls the driving force to the crystal phase and w controls the height of the double well. Here we have set w = 1.
In the specific case of a liquid-solid interphase, the parameter e(r) represents surface tension of the crystal that depends on the direction of the crystal orientation. In curvilinear coordinates, the crystal orientation is measured by the angle ψ that the crystal makes with thex 1 versor, i.e., tan(ψ) = where e 0 controls the magnitude of penalization in the free energy, n represents the anisotropy of the system (n = 4 for a square, n = 6 for a hexagon, and so on), and δ 0 controls the strength of the anisotropy. In this approach the dynamics of the phase transition can be studied through a relaxational equation of the form (at constant T): where µ φ is the mobility of the system and δF δφ is the functional derivative of F in terms of the order parameter φ. Defining the free energy F as F = L(φ, ▽φ, T)d 2 r, the functional derivative is written: where the curvilinear gradient and divergence are ▽ = g αβ ∂ α and ▽· = 1 √ g ∂ α √ g. In the case of adiabatic crystal growth, the effects of temperature variation in the neighbor of the liquid/crystal interface are key, and an equation for the diffusion of heat is also needed. Note that local changes in the temperature modifies ]. An evolution equation for the temperature field can be obtained from the conservation of enthalpy: where LB represents the Laplace-Beltrami operator Here K is an adimensional parameter, proportional to the latent heat and inversely proportional to the heat capacity of the solid phase. For the sake of simplicity, we establish the same constant of diffusion for the solid and liquid phases.
Since the driving force of the interphase is the subcooling, the parameter m must be a function of m(T e − T). In that way, we can model m as: where K 1 and K 2 are parameters such that K 1 < 1 thus |m| < 1 2 . Thus, in this model the adiabatic evolution of the system is given by the coupled differential Equations (5) and (7). To study crystal growth on the sphere we write the equations in spherical coordinates, where points on the sphere are specified by r = (θ , ϕ) where x 1 = θ ∈ [0, π] is the polar angle measured from a fixed zenith direction, and x 2 = ϕ ∈ [0, 2π) is the azimuth angle, Figure 1 shows the mesh on the sphere due to the spherical parametrization. In these coordinates, the arc length ds is given by ds 2 ≡ |dr| 2 = a 2 dθ 2 + sin(θ ) 2 a 2 dϕ 2 , and the differential operators take the form: The evolution equations are numerically solved by a finite difference scheme, forward in time and centered in space. As the initial crystal seed we use a geodesic circle. Multiple growing nuclei can be modeled by seeding different circular (geodesic) seeds on different parts of the sphere. Typically, we set the parameters e 0 = 0.01, µ φ = 1 0.0003 , δ 0 = 0.01 . To model the double potential well we set K 1 = 0.9 and K 2 = 10. Due to the spherical coordinates, our computational grid is not uniform, see an scheme in Figure 1, and the cell size decreases in the polar angle direction. At the pole of the sphere the element size is zero, and hence, the pole is singular point. Therefore, in our analysis, we ignore the growth near the poles. In order to use equivalent grids, we vary the number of grid points for spherical substrates of different radius (up to n θ × n ϕ = 512 × 1024 for a sphere of radius 4). The angular step, are then, θ = π n θ and ϕ = 2π n ϕ .

RESULTS
In order to test the model, we first consider isothermal crystallization on spheres. Figure 2 shows snapshots at different times of isothermal growth of nuclei for different degrees of anisotropy. From top to bottom the different rows in this figure show the evolution of the system for n = 1, n = 3, n = 4, and n = 6, respectively. As clearly seen from the panels, on spheres, initial circular nuclei acquire symmetric shapes as a consequence of the symmetry in the surface tension, similarly to isothermal crystallization in planar geometries (Kashchiev, 2000;Kelton and Greer, 2010). We now consider the features of dendritic growth on spherical substrates. Figure 3 shows the shape of dendritic nuclei grown on spherical substrates of different radius a = 1.5, 2, 3, 4. To make this comparison we use the same circular seed of geodesic radius r = 0.25, and the panels show the crystallites at the same time t = 0.1.
From this figure it can be qualitatively observed that the nuclei look less ramified on substrates of smaller size (compare Figures 3A,D). This suggests that curvature may originate a reduction in the ramifications of nuclei. In order to study this in more detail we calculate the degree of Circularity C of nuclei under same conditions, but on spherical substrates of different radius: where A is the area of the nucleus, and L its perimeter. Note that the Circularity of a planar circular nucleus is C ≡ 1, while the Circularity of a fractal nuclei approaches zero C ≡ 0. Figure 4 shows the Circularity of nuclei grown on different spherical substrates. Note that for different values of K, the circularity C continuously increases for smaller substrates, showing that the curvature contributes to the reduction of ramifications and the propagation of more compact nuclei. Phase field models can be sensitive to the discretization (mesh-induced anisotropy). In order to study if these results are influenced by the meshes, we have run simulations of growing crystal of different orientations (Dobravec et al., 2020). The variations observed in the circularity plotted as error bars in Figure 4, and are of the order of symbol size showing the robustness of the results. Panels (A-D) correspond to crystallites grown on substrates of radius a = 1.5, 2, 3, 4, respectively. Note that for small substrates (higher curvature) the dendrites look less ramified. The color bar shows the amplitude of the order parameter φ. It is also interesting to study how the substrate's curvature influences the selection of tip radius and propagation velocity of dendrites. In order to obtain the tip radius, we measured from the simulations the geodesic curvature of the interface at the tips as a function of the size of the spherical substrate. Suppose that the curvature of a curve C belonging to a surface S is k. The curvature of the curve at that point can be composed in two components, the normal curvature k n and the geodesic curvature k g (see Figure 5) (O'Neill, 1997). While the normal curvature is obtained by projecting curvature vector on the surface normal, the geodesic curvature is obtained by projecting on the tangential plane at that point. The curvatures are related by the identity: Note that the geodesic curvature is the intrinsic curvature of a curve inside the surface, and its the relevant quantity related to curvature-driven growth on surfaces. To measure the tip radii from the simulations, we first extracted the contour interface and calculated the curvature k at the tip by using the Frenet-Serret formulas. We then calculate the normal curvature k n by projecting to the sphere's normal and that point. Then, the geodesic curvature of the tip is obtained from Equation (13). Finally, the radius of curvature R g of the tip can be approximated FIGURE 5 | Sketch of the curvature of curve embedded in a surface in the 3D-space. This figure shows the curvature k on curve C which has a normal n on a surface S. The normal of the surface is N, and its projection with the curvature k is the normal curvature k n . The geodesic curvature k g , that is perpendicular to k n is, by definition, k 2 g = k 2 − k 2 n .
by k g ≈ 1 R g , which is valid for small circles on a sphere, as shown in Figure 6. Figure 7 shows how surface curvature affects the selection criterion. Here we plot the selected velocity of the tips as a function of the radius of the tip, for dendritic crystals growing on spheres of different curvature (the radius of the spherical substrates is indicated right on the side of the data). Note that as the size of the spherical substrate decreases, the selected tip radius increases. This is the same as reported above, the dendrites get more rounded for spherical substrates of higher curvature (smaller spheres). On the contrary, the selected velocity shows a non-monotonic behavior. For smaller spherical substrates the tip velocity increases at first, but later starts to decrease. It is interesting to note that a similar non-monotonic behavior of the tip velocity as a function of tip radius is also observed in planar geometries when a different degree of undercooling FIGURE 7 | Tip velocity as a function of the tip radius. This figure shows the tip velocity as function of the tip radius for dendrites growing on spheres of different radii (labeled at the right of each symbol). Errors bar are of the order of symbol size and omitted for clarity.   Panels (A,B) correspond to the crystal interface of isotropic and six-fold symmetry, respectively, as seen in the θ − ϕ plane (t = 0.116). Panel (C) shows the crystal interface deviation form a geodesic circle ǫ(ϕ) at t = 0.056 (below), the structure factor of this interface configuration (upper left), and the exponential growth of the peak of the structure factor (upper right). Panel (D) shows the circularity of the growing nuclei for a spherical substrate of different radius a, as obtained by the numerical solution of the interface model (t = 0.0001, and α = β = γ = 1). are used. Such behavior is roughly described by the spherical approximation which considers the tip of the dendrite as a growing sphere (Langer, 1980). Our simulations show that not only the undercooling, but also the substrate's curvature selects both, tip's radius and velocity.
In order to rationalize how the underlying curvature can modify the shape of dendrites on curved substrates, here we consider the interfacial instabilities associated with dendritic growth. To do this we use an equation describing the temporal evolution of the interface curvature: where V(κ g ) = κ g + ακ 2 g − βκ 3 g , and s and κ g are the arc length and the geodesic (intrinsic) curvature of the interface. Originally, this equation was used as a geometrical model for interface evolution in a variety of different problems (Brower et al., 1984;Kessler et al., 1984), including dendritic crystallization. Although in the original model the key property of the interface was the total curvature κ, here we slightly modify the model by using the geodesic curvature κ g , which is the curvature of the interface as seen from the curved substrate. In general, the geodesic curvature of a curve is also known as the intrinsic curvature of a curve (O'Neill, 1997).
We now use this model to develop an instability analysis of initial circular nuclei growing on a spherical substrate. This can be done by considering normal perturbations to an initial circular nuclei of geodesic radius r, such that the equation of the interface can be written as − → x (ϕ, t) = [r(t) + ǫ(ϕ, t)], where ǫ is the perturbation and ϕ is the azimuthal angle. In general, the geodesic curvature of a closed curve on the sphere can be written as (Parisio et al., 2001): a 2 cos (r/a) sin (r/a) 2 + 2 cos (r/a)( ∂r ∂ϕ ) 2 − a sin (r/a) ∂ 2 r ∂ϕ 2 [a 2 sin (r/a) 2 + ( ∂r ∂ϕ ) 2 ] 3/2 (15) where r is the local geodesic radius of the curve (the coordinates r, ϕ are usually known as geodesic polar coordinates; Gómez et al., 2015).
Thus, to first order in the perturbation ǫ, we can approximate: Now, by expanding the perturbation in the polar angle using the eigenfunctions ǫ m = ǫ 0 e λ m t cos (mϕ), we obtain the growth rate: We note that for large spherical substrates (r << a) the above expression reduces to growth rate obtained previously for dendritic growth in the plane λ 0 m = (m 2 − 1)[V ′ (1/r) − γ m 2 /r 2 ]/r 2 (Brower et al., 1984;Kessler et al., 1984).
Thus, the first correction of the amplification factor due to curvature can be written as λ m ≈ λ 0 m /(1 − r 2 /3a 2 ), such that although the growing nuclei are still linearly unstable on spheres, the exponential amplification of perturbations reduces for higher curvatures.
In order to study the evolution of the interface at early times (linear behavior), but also at long times (non-linear behavior), we numerically solve Equation (14), on spherical substrates of different sizes, by using a predictor-corrector approach. Figure 8A shows a typical interface morphology obtained at early times, as seen in the θ − ϕ space (t = 0.116). In Figure 8B we show that anisotropic dendritic growth on spheres can also be modeled by including a factor 1+ε cos(nθ ), multiplying the right hand of equation 14 (this figure corresponds to a simulation with hexagonal symmetry n = 6). In the lower panel of Figure 8C, we show the interface deviation from a geodesic circle ǫ, as a function of ϕ (t = 0.056). This can be used to analyze the evolution of the different modes during dendritic growth. In the top panels of Figure 8C we show the structure factor of the lower figure, and the evolution of the peak of the structure factor. As predicted by the linear analysis there is a range of unstable modes, and the amplification factor grows exponentially in time, until the non-linearities of the equation slows-down the growth, similarly as occur in other unstable systems like spinodal decomposition (Vega and Gómez, 2009).
It is worth to study the degree of ramifications as obtained during dendritic growth for this interface equation model. In Figure 8D, we show the circularity of growing nuclei, for spherical substrates of different sizes. Similarly to obtained with the phase field model, here we also observe a reduction of ramifications as the size of the substrate is reduced. Here, we note that a similar result was obtained for a phase field crystal model, where the degree of ramifications of growing nuclei was also observed to be reduced for smaller spheres. In such cases it was observed that growing nuclei evolve from a sixfold branching geometry to a ribbon as the size of the spherical surface is reduced (Köhler et al., 2016). Note that our results show that the reduction of ramification by curvature is also found even when the elastic interaction terms due to the crystal lattice are not taken into account, as in our model.
A similar mode expansion analysis was developed to study the instability of an interface separating two immiscible viscous fluids on a sphere (the Saffman-Taylor problem on a sphere) (Parisio et al., 2001). In such a work, a perturbation analysis up to second order in ǫ, showed an asymmetry on the tipsplitting behavior depending on the region where the interface was evolving. While tip-splitting is still present for interfaces on the northern hemisphere, tip-splitting completely vanishes in the southern hemisphere, where is completely replaced by a finger tip-sharpening mechanism.
Here we also found that tip-splitting is suppressed for large nuclei whose interface goes beyond the equator. Figure 9 shows the temporal evolution of a large crystal (in red). Note that the interface still becomes linearly unstable, but there is no tip-splitting leading to dendritic patterns. As can be observed from panel Figures 9B,C there is a tip-narrowing mechanism rather than tip splitting, similarly as predicted for fluids (Parisio et al., 2001). This is a simple direct consequence of the geometry and topology of the sphere, where such a large growing nuclei produce a converging flow, where the mechanism of tip-splitting is completely suppressed (Thomé et al., 1989).

CONCLUSIONS
In this work, we have explored the features of dendritic crystal growth on spherical substrates. In general, the underlying curvature of the substrates can be used to control the shape of crystals and propagation velocity. We found a reduction in the ramification of growing nuclei as a consequence of the curvature of spheres, which can be understood from a linear instability analysis of the crystal-liquid interface. For large nuclei whose size surpass half the surface of the sphere (the interface goes beyond the equator), there is a vanishing tip-splitting mechanism. For such nuclei the interface is still unstable, but the non-linear process of tip-splitting is replaced by a tip-narrowing behavior. Our results show that geometry and topology can drastically affect crystallization and aggregation mechanisms on spherical substrates. This could help to control the features of crystalline textures on micro/nano particles.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary materials, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
LG designed the research. LO performed the numerical simulations. Both authors discussed and analyzed the results and contributed to the manuscript preparation.

FUNDING
This work was supported by the National Research Council of Argentina, CONICET, the Fondo para la Investigación Cientifica y Tecnologica (FONCYT, Grant PICT-2017-3611), and Universidad Nacional del Sur.