Dendritic growth of ice crystals: a test of theory with experiments

Motivated by an important application of dendritic crystals in the form of an elliptical paraboloid, which widely spread in nature (ice crystals), we develop here the selection theory of their stable growth mode. This theory enables us to separately define the tip velocity of dendrites and their tip diameter as functions of the melt undercooling. This, in turn, makes it possible to judge the microstructure of the material obtained as a result of the crystallization process. So, in the first instance, the steady-state analytical solution that describes the growth of such dendrites in undercooled one-component liquids is found. Then a system of equations consisting of the selection criterion and the undercooling balance that describes a stable growth mode of elliptical dendrites is formulated and analyzed. Three parametric solutions of this system are deduced in an explicit form. Our calculations based on these solutions demonstrate that the theoretical predictions are in good agreement with experimental data for ice dendrites growing at small undercoolings in pure water.


Introduction
It is well-known that the processes of phase transformation from the metastable and nonequilibrium states underlie the preparation of many materials with specific physicochemical properties. As this takes place, the growth features of dendritic structures evolving in undercooled liquids determine the emerging microstructure of materials obtained in solidification processes [1][2][3][4][5][6]. For a theoretical description of the stable growth mode of dendritic crystals, the theory of microscopic * Author to whom any correspondence should be addressed. solvability was originally developed for pure (one-component) melts [7][8][9], which makes it possible to select the relationship between the growth velocity of a dendrite tip, its radius of curvature, and melt undercooling. This theory is based on two nonlinear equations-the selection criterion and the undercooling balance. Using these algebraic equations, it is possible to independently determine the velocity and radius of curvature of the dendritic tip depending on the system undercooling (the driving force of dendritic growth). Then this theory was extended to binary melts [10,11], dendritic growth under conditions of convective fluid flow [12][13][14][15][16][17], and also to processes of local nonequilibrium (rapid) crystallization [18][19][20]. 3D phase-field simulation of non-axisymmetric ice dendrite growth [21,22]. From left to right: stellar dendrite with secondary branching, stellar dendrite, fern-like dendrite, and 12-arms star. Reproduced from [22]. CC BY 4.0.
Aside from [21,22] (see figure 1), the vast majority of theoretical and numerical studies using solvability theory are based on the assumption of axisymmetric dendrite tip (for instance [23]). However, experimental data show that dendritic crystals often represent non-axisymmetric structures, i.e. their tip shapes may have different curvature radii in the basal and perpendicular planes (ice crystals represent a good example of such non-axisymmetric dendrites) [24][25][26][27][28]. These experimental observations have yet never been ensured by a theory that allows selecting a stable growth mode for such crystals. The present work is devoted to the selection theory of nonaxisymmetric dendritic crystals growing with various crystalline symmetries in the basal and perpendicular planes. Our approach is based on the analysis of two-dimensional crystal growth, taking into account the fact that the dendritic tip in these planes evolves with the same velocity. The developed theory is tested against experimental data that confirm our theoretical conclusions.

Analytical theory
Let us consider the stationary growth of dendrite in the form of an elliptical paraboloid from an undercooled one-component melt (see figure 2). For the sake of simplicity, we neglect the melt flow that could be induced by the temperature and pressure gradients. The elliptical paraboloid is defined by two forming parabolas in perpendicular planes with diameters ρ 6 and ρ 2 . Here we mean the growth of ice crystals having sixfold symmetry in the basal plane and two-fold symmetry in the perpendicular plane. Keeping this in mind we use subscripts 6 and 2. Noteworthy, the present theory can be generalized to the growth of a three-dimensional dendritic crystal with another crystalline symmetry (for instance, 6 and 3). Next, assuming that the dendritic interface is isothermal, we describe the temperature field in the melt by a single coordinate that specifies the isothermal surfaces and represents an elliptical paraboloid [29].
In this case, the interface equation in dimensionless coordinates x, y, z takes the form where the length scale is a 2D T /V, D T is the thermal diffusivity, V is a constant growth rate, ω is an isothermal level set variable and b = (ρ 6 − ρ 2 )D T /V is a dimensionless parameter of ellipticity. Note that in the case of stationary growth, the dendritic shape is fixed, and the ratio k of its tip diameters is constant,

Selection criterion
The selection criterion allows us to obtain the dendritic growth velocity and tip diameters in two planes of the hexagonalpacked crystal lattice, having the six-fold symmetry in its basal plane and two-fold symmetry in the perpendicular vertical plane. In the 2D and axisymmetric 3D cases, this criterion has been found both analytically and numerically [7,11]. For the nonaxisymmetric 3D problem, the axisymmetric solution can be used as an approximation in the vicinity of dendritic tip [30]. Thus, to find three parameters, namely the growth velocity V and two dendrite tip diameters ρ n (n = 2, 6), we use two solvability criteria in the two perpendicular planes of the dendrite (basal, where n = 6 and vertical, where n = 2). Therefore, let us write out the general form of 2D selection criterion for a single-component system without convection as [11][12][13][14][18][19][20] where d 0 is the capillary length constant, σ 0n is the solvability constant, α dn is the anisotropy parameter (stiffness), n is the order of symmetry (it indicates a concrete crystalline plane), P gn = (ρ n V)/(2D T ) is the symmetry-dependent Péclet number and a 1n = 8σ 0n 7 3 56 (see, for more details, references [13,26,31]). It is significant to note that the selection criterion (3) connects the averaged tip diameter ρ = ρ 6 − b = ρ 2 + b and the tip velocity (the rhs also depends on these unknowns through P gn ).

Undercooling balance
The second relation between the growth velocity and tip radii can be retrieved from the undercooling balance representing the driving force of crystal growth. The total undercooling ΔT = T m − T ∞ at the dendrite tip consists in several contributions where T m is the melting temperature of a planar front, T ∞ is the temperature far from the growing dendrite, and ΔT T stands for the thermal undercooling, which reads as (see appendix) Here T Q = Q/c p represents the adiabatic temperature, Q is the latent heat of solidification, c p is the specific heat, P T = (ρ 6 + ρ 2 )V/(4D T ) is the Péclet number based on the averaged tip diameter ρ = (ρ 6 + ρ 2 )/2. The curvature contribution ΔT R stands for the undercooling induced by the Gibbs-Thomson effect, which can be written in the form of Note that the total undercooling balance (4) connects ρ = ρ 6 − b = ρ 2 + b, V and ΔT.
Using the stability criterion in two perpendicular planes with a different order of symmetry together with the undercooling balance allows us to have a closed system of three equations with three unknowns (two radii ρ 6 and ρ 2 of forming parabolas and constant growth velocity V).

Exact analytical solutions
Let us find now the exact analytical solutions of the nonlinear set of three governing equations (3) and (4) in a parametric form. Note that equation (3) represents two selection criteria for the larger (ρ 6 ) and smaller (ρ 2 ) dendritic tip diameters (ρ 6 > ρ 2 ). To do this, we first equate the growth rate V expressed in terms of ρ 6 and ρ 2 (see criterion (3)) to each other and obtain V(ρ 2 ) and ρ 6 (ρ 2 ). Finally, substituting them into the undercooling balance (4), we get the first analytical solution where P T (ρ 2 ) = (ρ 6 (ρ 2 ) + ρ 2 )V(ρ 2 )/(4D T ), parameter n equals to 6 or 2, and ρ 2 plays the role of a decision variable (free parameter that can be tunned). Expressions (7) represent the first parametric solution.

Experimental tests
Let us compare the present theory (expressions (8)) with experimental data [24,25] on the growth of ice crystals (physical parameters are given in table 1). Figures 3-5 illustrate the dendritic tip velocity V and its tip radii in the basal (larger radius ρ 6 ) and perpendicular to it (smaller radius ρ 2 ) planes. As is easily seen, the tip velocity monotonously increases and tip radii monotonously decrease when increasing the melt   undercooling ΔT. The theory and experimental data are in good agreement in the range of small and moderate undercooling.
The smaller dendrite tip diameter ρ 2 was not directly measured but calculated following the results by Furukawa and Shimada [25,32] as the direct ρ 2 data could not be obtained in the ISS (International Space Station) experiments. So, their semi-empirical law takes the form where Δ = ΔT/T Q is the dimensionless undercooling. An important point is that this formula is a generalization of experimental data. The smaller tip diameter ρ 2 calculated on the basis of our theory (expressions (8)) as a function of the melt undercooling is compared with experimental data (expression (9)) in figure 5. A slight deviation of the curves is explained by the fact that the attitude of the dendrite could not be easily controlled in experiments [25].

Conclusion
In summary, the theory of stable three-dimensional dendritic growth in the form of an elliptical paraboloid is developed on the basis of undercooling balance condition and two stability criteria written out in two perpendicular dendritic planes (the basal plane described by the larger tip diameter ρ 6 and perpendicular to it plane described by the smaller tip diameter ρ 2 ). This theory leads to three non-linear algebraic equations determining three growth parameters-the dendrite tip velocity V and its tip diameters ρ 6 and ρ 2 as functions of the melt undercooling. The solution of these equations governing the stable growth mode is analytically found and tested against experimental data for ice crystals. We demonstrate that predictions of our theory agree well with experimental data in the range of small undercoolings. The theory is applicable for dendrites with non-axisymmetric morphology such as ice crystals with tip shapes of elliptical paraboloids.

Acknowledgments
The present research work consists of computational and theoretical parts, which were supported by different financial sources. In this connection, DVA acknowledges support

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

Appendix. Boundary integrals for elliptical dendrites
Let us introduce how to derive the Horvay and Cahn solution describing the elliptical paraboloid from the boundary integral method. Horvay and Cahn [29] obtained the shape of an isothermal dendrite having an elliptical cross section. In this case the surface of the dendrite in dimensionless coordinates has the form (1). The dimensionless undercooling as a function of P T can be written as [33] Δ − d c ρ where (in dimensionless variables) Here Δ = (T m − T ∞ )/T Q is the dimensionless undercooling, β is the kinetic coefficient, variables x, x 1 and ζ are measured in units of ρ/P T , and variables t and τ are measured in units of ρ/(P T V). Further, the substitution of interfacial function ζ = z from (1) leads to Replacing two integration variables τ and y 1 by ω 1 and z 1 respectively as we can obtain Integrating the right-hand side of equation (14) over x 1 and z 1 yields Where it has been taken into account the following [34] ∞ 0 exp(−μ 2 u 2 )du u 2 + β 2 = π 2β erfc(βμ) exp(β 2 μ 2 ).
Replacing the variable ω 1 by u = 1 + ω 1 ω−b , we obtain where Next, applying the method of differentiation and integration J(ω) and considering that we arrive at Keeping in mind that at the interface ω = P T we obtain For b = 0, the cross-section is circular and one can obtain the Ivantsov's solution [13,35].