Dynamical scaling for underdamped strain order parameters
quenched below first-order phase transitions

Interacting systems such as magnets, binary fluids, liquid crystals, and quantum spin models can be probed by the dynamical evolutions of order parameters, after temperature or coupling-constant quenches below transition [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. Such phase ordering can be quantitatively described by a time-dependent, two-point, order-parameter correlation $C(\vec{R},t)$, that can exhibit dynamical scaling, dependent on space and time through a single scaled variable[1, 5, 6, 7, 8, 9, 10, 11, 12, 13] ${\bar{R}}\equiv|\vec{R}|/L(t)$, with consequent data collapse onto a single scaled curve $G(\bar{R})$. The coarsening length $L(t)$ is a measure of the typical spacing between domain walls separating competing order parameter (OP) phases. It can increase as $L(t)\sim t^{\alpha}$, where the exponent $\alpha$ is independent of material parameter values, but could depend on the nature of the OP dynamics, the number of components of the order parameter $N_{\rm OP}$, the number of competing low-temperature variants $N_{V}$, and the spatial dimension $d$. There can be a sequential appearance of different exponents, during coarsening [1].

In various dynamic models, the exponents $\alpha$ have been estimated by heuristic arguments, energy dissipation matchings, Gaussian fluctuations of domain wall profiles, self-consistent correlation function dynamics, and through numerical simulations [1, 5, 6, 7, 8, 9, 10, 11, 12]. The Allen-Cahn relaxation equation [1] for a nonconserved OP and the Cahn-Hilliard equation [1, 4] for a locally conserved OP are familiar models, both with a single time derivative, and typically use a double-well Landau free energy describing a second-order transition. It is of much interest to explore other phase-ordering dynamics; and to develop systematic methods of estimating coarsening-length exponents.

Solid-solid structural transitions have $N_{\rm OP}$ strain-tensor components as the order parameters [14, 15, 16, 17, 18, 19, 20, 21], with Landau free energies having $N_{V}$ competing minima, for the $N_{V}$ different unit cells, or ’variants’. The high-temperature, high-symmetry crystal structure is “austenite,” and the low-temperature, low-symmetry structures are “martensite.” Martensitic transitions can be described by a Bales-Gooding or BG-type strain dynamics [19] that has several features [19, 20], which differ from the more familiar magnetism-inspired phase orderings. Firstly, the Landau free energy [16] can have triple wells in the OP, describing a first-order phase transition, with a minimum also at zero values of the OP. Secondly, the dynamics is underdamped, with a Newtonian inertial term or double time derivative, describing acceleration of the order parameter. Thirdly, there is global momentum conservation, with the single time derivative damping term, suppressed at long wavelengths[19, 20]. Fourthly, with a 2D dynamics [20] generalized from [19] 1D, the order parameters have an additional power-law anisotropic interaction, coming from an elastic St-Venant compatibility constraint [15, 16], as used in several contexts [21].

Monte Carlo simulations in 2D of a discretized strain-pseudospin martensitic model Hamiltonian with power-law anisotropic (PLA) interactions show interesting evolutions under a temperature quench. For example, for successive quenches approaching the transition from below, the conversion time from seeded austenite evolving to martensite domains rises sharply, with these time delays caused by entropy barriers [18]. Clearly, martensites with continuous-strain dynamics [19, 20] are worth examining, in the framework of phase ordering ideas [1]. We here focus on effects of the triple well Landau term in the underdamped dynamics, and suppress the power-law anisotropic interaction, which will be considered in a subsequent publication.

In the first part of this paper, we apply phase-ordering ideas to the BG strain dynamics, with $d=2$, $N_{\rm OP}=1$, and $N_{V}=2$, in a sixth-order Landau polynomial in the scalar OP strain. We consider only Landau and Ginzburg terms in the OP dynamics to numerically determine the dynamic structure factor or OP-OP correlation, finding dynamical scaling in a coarsening length $L(t)\sim t^{\alpha}$. The exponent takes on sequential values such as $\alpha=2/3,1/2$ over time windows, whose widths depend on the quench temperature $T$. For deep quenches, there is an exponent $\alpha=1/3$, and a final $\alpha\simeq 0$ flattening to a constant, analogous to ’coarsening arrest’ [22]. These 2D results are found to persist, for 3D.

In the second part of the paper, we use a scaled form of the underdamped OP dynamics, to obtain a dynamics for the OP-OP correlation $C(R,t)$. Inserting the dynamic scaling form, $C=G(R/L(t))$, yields a nonlinear, underdamped kinetics for the curvature or inverse coarsening length $g(t)\equiv 1/L(t)$, with coefficients evaluated at dominant coarsening-front separation. Here, to close what would otherwise be an infinite hierarchy, a spatial average of internal domain-wall factors is made, reducing the correlation between the chemical potential and order parameter, to the OP-OP correlation $G(gR)$. The curvature kinetics solutions from balancing kinematic and force terms are simple power-law decays $g(t)\sim 1/t^{\alpha}$ in sequential time windows, showing exponents $\alpha=2/3,1/2$ values, independent of $d$, in agreement with simulations. For deep quenches, a toy model including higher powers of the curvature, explains the $\alpha=0$ coarsening arrest, as a metastable trapping of curvature to a nonzero value.

As a check, the curvature kinetics method is applied to Cahn-Hilliard dynamics, yielding the well-known [1, 8, 9] values of $\alpha=1/3$ for a scalar order parameter $N_{\rm OP}=1$, and $\alpha=1/4$ for vector order parameters with $N_{\rm OP}\geq 2$, all independent of $d$.

The plan of the paper is as follows. In Sec. II, we state the Cahn-Hilliard and Bales-Gooding types of OP dynamics, and scale to absorb all, or most, of the OP temperature dependencies. Section III defines the OP-OP correlation functions and their dynamical scaling. Section IV shows the coarsening textures, numerically demonstrates dynamical scaling of the two-point correlation function, and states the obtained exponents. The second part of the paper, starting in Sec. V, obtains the correlation function dynamics, and inserts a dynamical-scaling ansatz. Section VI extracts the curvature kinetics, and predicts power-law exponents $\alpha$ that match the numerics. Finally, Sec. VII contains a discussion and comments on future work. Details are given, in Appendix A, of the closure approximation; in Appendix B, of coefficient signs in the curvature kinetics; in Appendix C, of estimation of coarsening exponents.

In this section, we define the OP-OP and related correlations, and their dynamical scaling forms.

For a one-component martensitic-strain order parameter $e(\vec{r},t)$, we consider a two-point correlation between OP’s at $\vec{r}=\vec{R}+\vec{r}_{0}$, and $\vec{r^{\prime}}=\vec{r}_{0}$ on a $d$-dimensional lattice, and at equal times $t$. With an average over all origins ${\vec{r}_{0}}$ and over many runs, OP-OP correlations are dependent only on the separation $\vec{R}=\vec{r}-\vec{r^{\prime}}$,

With a Fourier expansion $e(\vec{r},t)=\frac{1}{\sqrt{N}}\sum_{\vec{k}}e^{i\vec{k}\cdot\vec{r}}e(\vec{k},t)$, we get

where the time-dependent structure factor is

Since the OP is real, its Fourier coefficients $e(\vec{k},t)^{*}=e(-\vec{k},t)$ and so $S(\vec{k},t)=S(-\vec{k},t)$.

Another correlation that enters is the chemical potential-order parameter or $\mu$-OP correlation:

With a Fourier expansion $\mu(\vec{r},t)=\frac{1}{\sqrt{N}}\sum_{\vec{k}}e^{i\vec{k}\cdot\vec{r}}\mu(\vec{k},t)$, we get

where only the symmetric part survives, and the $\mu$-OP Fourier correlation is the real part,

Since the anisotropic compatibility interaction is switched off, the system is isotropic, and we can work throughout with averages $<...>$ that now include angular averages, so correlations become $C(\vec{R},t)\rightarrow C(R,t)$ and $S(\vec{k},t)\rightarrow S(k,t)$,$S^{(\mu)}(\vec{k},t)\rightarrow S^{(\mu)}(k,t)$ where $R\equiv|\vec{R}|$, and $k\equiv|\vec{k}|$.

The OP-OP correlation $C(R,t)$ for a given time $t$ will show a fall-off in separation $R$; while at a given separation $R$, it will increase with time. Since $\sim L(t)$ is the scale of the OP-correlation region, it will also increase with $t$. As $L(t)$ is also the separation of the OP-dips at domain walls, these strain patterns must coarsen with time.

Dynamical scaling says that [1, 13] (i) the time $t$ enters only through the length $L(t)$; (ii) the separation $R$ appears only in scaled form as $\bar{R}\equiv R/L(t)$, [and in Fourier space, the wave vector $k$ appears only as $\bar{k}\equiv kL(t)$]. Further, in a common assumption, (iii) the coarsening length scale grows as a power law $L\sim t^{\alpha}$, which is plausible [1], but here is justified. The OP-OP correlation is then

We find later that the zero-separation values, in a few hundred time steps, on the onset of dynamical scaling, are insensitive to time, $C(0,t)\simeq C(0,t_{\rm onset})$, so we henceforth suppress the constant denominator.

In Fourier space, the scaling behaviour is

For sharp domain walls, of widths small compared to $L(t)$, Porod’s law holds for the Fourier space structure factor [1, 7]

For $N_{\rm OP}=1$, we have $\chi\simeq 1/(kL)^{d+1}$, and a log-log plot of $S(k,t)$ versus $k$ will be a straight line with slope $-(d+1)$, with the length $L(t)$ extracted from the intercepts $\ln\chi=-(d+1)\ln k-\ln L(t)$. Then with this scaling length, replots of $\ln\chi$ versus $\ln(\bar{k})$, as well as $G(\bar{R})$ versus $\bar{R}$, will show data collapse of different-time curves.

The coarsening curvature is defined [1] as $(d-1)/L(t)$, but we will simply work with an inverse length, or

and call this the “curvature,” as it indeed is, for $d=2$. Thus the scaled variables are

We later show that the correlation dynamics results in a nonlinear kinetics for the curvature $g(t)$, that has power-law solutions $g\sim 1/t^{\alpha}$, explaining the observed coarsening-length behaviour $L\sim t^{\alpha}$.

In this section, we present numerical work, showing domain-wall (DW) textures, demonstrating dynamical scaling, and extracting exponent behaviour of coarsening curvatures. In an actual experimental quench, the physical temperature is changed suddenly, and the coarsening probe is followed in physical time. Any universal curve independent of quench temperatures $T$, is obtained only by later scaling the measured physical values, by factors dependent on $T$. As mentioned, we mimic this experimental procedure in simulations by using the unscaled forms of both the CH and BG dynamics in (2.3)and (2.10); only then doing scaling on the numerical data obtained.

The initial high-temperature state for both the CH and BG cases, corresponds to a single well, with a minimum at $e(\vec{r},t=0)=0$, plus a few local fluctuations of the ordered variants. A numerical “quench” then corresponds to evolving at some suddenly lower temperature $T$. There will be an early-time regime where the ordered phases expand, to crowd out the initially dominant $e=0$ background as a DW vapour phase; while at intermediate and long times, there will be only competing variants, separated by domain walls, as a DW liquid phase.

We thus start with a dilute set of initial seeds of both martensite variants equally, in a sea of $e=0$ austenite. The martensite conversion fraction $n_{m}(t)\equiv\frac{1}{N}\sum_{\vec{r}}e^{2}(\vec{r},t)/|\bar{\varepsilon}|$ is initially nonzero only at $2\%$ random sites surrounded by $2\times 2$ unit cells where $e=\pm{\bar{\varepsilon}}$, with $e=0$ elsewhere, and so $n_{m}(0)=0.08$.

We use an Euler-discretized dynamics of (2.3), (2.10) to find the OP evolution of $e(\vec{r},t)$, focusing on $\Lambda=1$. We take time steps $\Delta t=0.05$; and spatial derivatives $\nabla_{\nu}$ as finite-difference operators $(a_{0})^{-1}\Delta_{\nu}$ on a square unit lattice. For wave vectors $\vec{k}$ in the Brillouin zone, $\Delta_{\nu}\rightarrow 2\sin(k_{\nu}a_{0}/2)$ , with the grid scale set to $a_{0}=1$. A fast Fourier transform yields the Fourier coefficient $e(\vec{k},t)$; and hence the angularly averaged dynamic structure factor $S(k,t)=<|e(\vec{k},t)|^{2}>$ of (3.2). (To focus on DW time scales, a “strain-hardening” procedure is used [9].) Run averages are taken, over $N_{run}=5$. A reverse FFT yields the OP-OP correlations $C(R,t)$ of (3.1). The quench temperature $T$ is held fixed, for a holding time $t_{h}=20,000$ steps. The 2D system is of size $L_{\rm sys}^{2}=8192^{2}$.

Figure 2(a) shows the BG case single-run martensite or ordered-fraction $n_{m}(t)$ versus scaled time $t\bar{\varepsilon}^{4}(\tau)$. For temperatures just below $T=T_{4}$ or $\tau=\tau(T_{4})=0.74$ the martensite fraction rises slowly towards unity, with larger delays, closer to $T_{4}$. Above $T_{4}$ there is no conversion at all, as the initial martensite seeds dissolve back into the four-fold symmetry austenite. We will quench to temperatures sufficiently below $T_{4}$ so conversion delays are small.

Figure 2(b) shows the CH case ordered-fraction $n_{m}(t)$ versus scaled time $t/t_{\rm sc}=t\epsilon^{2}$, where its slow rise for $T$ close to $T_{c}$ is due to critical slowing down near a second-order transition. As will be seen, these intermediate times can nonetheless show exponent behaviour.

Figure 3 shows snapshots of evolving relief plots of the OP $e(\vec{r},t)$ versus position $(x,y)$, for early time evolutions for $t=100,250,1500$, and deep quenches. The strain textures or DW patterns, clearly coarsen with time.

Figure 4 again shows snapshots of coarsening, but now as evolving contour plots, for shallow, moderate and deep quenches. We will later consider several such quenches just below $\tau=+0.6,+0.3,+0.05$, corresponding to $\Delta\tau(T)\equiv\tau(T)-\tau(T_{4})=-0.1,-0.5,-0.8$. The first row is for a shallow quench, to $\tau=0.6$, and shows, in a background of $e=0$ austenite, many whorl-texture droplets at early times, like a DW vapour. The second and third rows are for moderate and deep quenches around $\tau=0.2,-0.1$, with wandering martensite-martensite interfaces, like a DW liquid. For deep quenches, the coarsening seems to slow, or be arrested to form a DW Glass, as discussed later.

As a benchmark, we start with the familiar CH case, in 2D and with a scalar $N_{OP}=1$ order parameter. Fig 5 shows the well-known results of dynamical scaling.

Figure 5(a) shows $C(R,t)$ versus $R$ curves at different times for a given quench, with $\epsilon(T)=-0.5$ shown. Extracting the coarsening length and re-plotting, we find (i) data collapse of different- time curves on to a common scaling curve $C(R,t)=G(g(t)R)$, for that quench; further, (ii) data collapse of different-quench scaling curves $G(\bar{R})$ onto a single $T$-independent curve. This is consistent with (2.14), which predicted a OP-scaled CH dynamics would be independent of temperature. The curvature exponents are consistent with the literature, with a long-time exponent [1] of $\alpha=1/3$ and, for temperatures close to $T_{c}$, an intermediate-time exponent [5] $\alpha=1/4$. Simulations in $d=3$, show [8] the same asymptotic $\alpha=1/3$ exponent, that is, thus independent of spatial dimension.

Figure 6 shows CH case log–log plots of the scaled curvature $g/{\bar{\varepsilon}}=g/|\epsilon|^{1/2}$ versus scaled time $t{\bar{\varepsilon}}^{4}=t|\epsilon|^{2}$, showing data collapse, not only for all times, but also for all temperatures, as expected from the OP-scaled form of (2.14). Existing results [1, 8, 9] have temperature-independent exponents at intermediate-times where $g(t)\xi_{\rm GL}(T)>1$ as $\alpha=1/4$; and a long times where $g\xi_{\rm GL}(T)<1$ as the Lifshitz-Slyozov exponent of $\alpha=1/3$. Extracting exponent values from data (as for BG case of Appendix C), we find averaged exponents and standard deviations as

close to $1/4$ and $1/3$, as in previous work.

Turning to the BG case quenches, Bales and Gooding in 1D find a backward-bending phase boundary of inverse damping versus temperature. Figure 7 shows that in the 2D case, the phase diagram of inverse-damping $\sqrt{\Lambda}\sim 1/\gamma$ versus temperature $\tau(T)$ also has the transition occurring at lower $\tau$, for decreasing damping. Below the boundary, the dilute initial seeds with small $n_{m}$, evolve to $n_{m}\rightarrow 1$ (martensite), while above it, $n_{m}\rightarrow 0$ (austenite). All quenches are to below the phase boundary.

For $d=2$ in previous Monte Carlo simulations of a related model, a complex textural energy was parametrized by a surrogate-droplet energy that was a universal inverted parabola versus a scaled evolving textural parameter [18]. As a check, our BG case dynamical evolutions were benchmarked against this energy parametrization [23].

The BG case plots of Figs. 8-10 show numerical results analogous to the CH case, demonstrating dynamical scaling, for the three quench regimes each, just below $\tau=0.6,0.3,0.05$, as mentioned. More in detail, all the Figures show the following: (a) correlations $C(R,t)$ versus $R$ curves at different times, for quench to a given $\tau(T)$; (b) scaled structure factor ${\chi}(kL)$, using the extracted $L(t)$, versus $\bar{k}=kL$ showing Porod’s tail behaviour of exponent -3 and Fourier data collapse; (c) dynamical scaling with data collapse of different-time curves on to a common scaling curve $C(R,t)=G(g(t)R)$, for that $\tau$ quench. However, the scaling curves $G(\bar{R})$, are different for different $\tau$, especially for the shallow quenches of Fig. 8; while for deeper quenches of Figs. 9 and 10, the curves are closer. This is consistent with the OP-scaled result of (2.15), that shows a residual $T$-dependence of the BG dynamics, through $\eta_{sc}(T)=\tau(T)/3{\bar{\varepsilon}}(T)^{4}$, that is, insensitive to $\tau$ at low temperatures. Also they show (d) log-log plots of the coarsening curvature versus time, showing different indicated exponents in different time windows, within the holding time $t_{h}$.

Figure 8 shows $\alpha=2/3$; Fig. 9 shows $\alpha=2/3$ followed by $1/2$; Fig. 10 shows $\alpha=1/2$ followed by $\alpha=1/3$, and a peculiar flattening to ‘$\alpha=0$’ at low temperatures. The trapped curvature $g_{0}=0.1$ is not just a finite-size effect, as for our system sizes, $g_{0}\gg 1/L_{\rm sys}\sim 10^{-4}$ .

The exponent mean values and standard deviations, with simple arithmetic average over all temperatures, are

that are close to $\alpha=2/3$ and $\alpha=1/2$. See Appendix C for a time-window procedure for extracting exponents.

All this is for 2D. We have also considered 3D coarsening textures, as shown in surface contour plots of Fig. 11. The 3D case also shows dynamical scaling as shown in Fig. 12, with the Porod’s law exponent now $-(d+1)=-4$, and data collapse as before. For successively deeper quenches, the exponents are found to be again close to $\alpha=2/3,1/2,1/3,0$, so the coarsening exponents are independent of spatial dimension.

The exponent behaviour has been extracted from various OP dynamics by many authors [1, 5, 6, 11], including Siggia through insightful heuristic arguments; Bray and Rutenberg through matching global and local dissipation; Ohta, Jasnow, and Kawasaki; Bray, and Puri, and Mazenko, through a fluctuating domain wall approach; and Langer et al. through a self-consistent truncation of correlation function dynamics. We suggest a complementary, and systematic curvature kinetics approach, to extract exponent behaviour.

In our theoretical analysis, we will use throughout, the OP-scaled form of the CH and BG dynamics, that fully or mostly, scales out the OP temperature dependence. We first obtain the dynamics of the two-point OP-OP correlation as done previously, elsewhere [5].