A mathematical model for assessing the impact of dual-level toxicity on aquatic biospecies and its optimal control analysis

1.   Introduction

The development of the thermodynamically based calculation methods consisting of the solution of the heat and mass transfer equations with the phase evolution equations has made significant progress in the creation of the production technologies for metallic alloys with specific properties ^[1,2]. Despite significant progress, the study and description of the fundamental physical processes that occur during solidification remains a relevant and important task of modern condensed matter physics ^[3,4,5,6]. Modern technologies of the processing of metals and their alloys allow us to achieve experimentally significant undercoolings, temperature, and concentration gradients, and, as a result, rapid boundary interfaces ^[1]. Under such specific conditions, a liquid can solidify not only into an equilibrium phase, but also into a metastable crystalline, or into an amorphous metallic glass state. Experimental results present, that such mixed states can significantly improve the physical properties of materials ^[7]. Difficulties of describing such effects are related to the lack of understanding of the physics of the glass transition ^[8,9], as well as to the difficulty of going beyond the equilibrium statistical physics. The latter is a consequence of the ergodicity violation of the amorphous melt due to the structural frustration ^[10,11].

In this paper, we have tried to overcome these difficulties using the gauge theory of glass transition ^[12], which considers topologically stable distortions in forming the partially ordered phase as structural units. It allows us to describe the liquid-glass transition using continuous field theory in the presence of the time-dependent temperature field. This approach will potentially benefit to the development of phase field theory for the glass-forming alloys.

The paper is organized as follows: (ⅰ) After the introduction in Section 1.1, a brief review of the existing approaches of modeling amorphization and vitrification in the phase field method is given; (ⅱ) the assumptions of the theoretical model considering the topologically protected perturbations as the main structural elements are formulated in Sections 2.1 and 2.2; (ⅲ) the system of differential equations, considering the heat transfer is derived, and its numerical implementation is described in Sections 3.1 and 3.2; (ⅳ) in Sections 3.3, we outline the details of the performed simulations and their relation to the real experiments; (ⅴ) results and discussion are summarized in Section 4. We discuss each set of simulations sequentially: Section 4.1—annealing from the disordered state at a constant temperature gradient for a certain time; Section 4.2—effect of changing of the shear modulus and variable thermal conductivity coefficient of heat transfer equation on the final structure, and correlation functions; Section 4.3—effect of the cooling rate fixed in time and varied along the $x$-axis and its qualitative influence on the final structure obtained after solidification; (ⅵ) conclusions in Section 5 outline the obtained results.

1.1.   Amorphization and vitrification in phase field models

The phase field (PF) and phase field crystal (PFC) methods ^[2] are widely introduced in materials science, in fact, they can be considered as a generalization of the mean-field approximation of the phase transitions theory ^[13,14,15]. The combination of these methods with nonequilibrium thermodynamics allows one to successfully describe rapid first order phase transitions ^[16,17]. The PF also links the nano- and mesoscopic levels for the description of the structural changes ^[18,19]. The fundamental challenges of the PF and PFC methods have either been resolved or, at least, are clear ^[20,21,22]. However, their application to the formation of non-ergodic amorphous phases has not been thoroughly investigated, since it requires going beyond the approximation of the mean field.

The application of such models to the non-ergodic amorphous phase states is not yet entirely clarified, since it requires a description of the separate ordering field corresponding to the glassy state, which in the general case must take into account the thermal conditions. The multi-scale combination of the thermodynamic approach with the microscopic dynamics of the particles (or defects) driven by the fluctuations is a known problem partially solved with the Smoluchowski equation ^[23,24]; however, it faces significant difficulties in accurate temperature accounting ^[25]. In particular, an attempt to describe the transition to a glass-like state using the PFC method was made in ^[26,27,28] for pure single component systems, and the manuscript presents the formation of fine-grained glass-like structure at high cooling rates. One can observe the sequential delayed crystallization (and recrystallization) from amorphous nuclei to the crystalline bulk in PFC ^[29,30]. In ^[31,32], a hybrid model consisting of a combination of a continuous mean-field model and a multiphase PF model was developed. There, the liquid-glass transition was considered as a noise-induced structural relaxation, taking into account the phase-interface slowdown.

The transition to a glass-like state can be carried out using the kinetic models for diffusion coefficients as made in ^[33,34,35] with the PF model for single-component and binary systems; in these works the glass-forming phase is considered as a bulk phase with a low diffusion strongly dependent on temperature. However, this model did not take into account the frustration, and non-ergodic effects inherent in glass, which are necessary to describe the formation of a fine-grained glass-like structure at high cooling rates. In a number of papers, disordered amorphous phases are considered as analogues of liquids, where crystallization (phase transition) is the same as the transition from liquid or amorphous phase and corresponds to the first order phase transition ^[36,37].

In ^[38,39,40], the glass transition was considered as a full-fledged phase transition, and the parameters of glassy state are considered in a bulk using the thermodynamic data obtained experimentally or with several approximations. One can supply the possessed challenges with some additional problems arising from the glass-liquid-solid interfaces description ^[32,41], confined amorphous states ^[42], and from the presence of the long-time relaxations leading to a nonlinear mobility of the glass-forming alloys ^[43].

The PF-models usually operate with the scalar order parameters, but the glass transition introduces the infinite landscape of the possible metastable states, which can be introduced by the uncertainty of the direction of the vector order parameter. In the present work, we do not pretend to the new phase-field interpretation, although there is a possible benefit in addition to the PF-method such as a continuous model of the kinetically sensitive glass formation model. We propose a new continuous model unifying the simultaneous description of kinetically sensitive formation of crystalline and amorphous phases which, due to its simplicity, can significantly improve the PF-method in the description of non-ergodic effects of glass formation.

2.   Gauge model of vitrification

2.1.   Assumptions and model analysis

In order to create a phase-field model of the solidification process, with the competition between crystallization and vitrification, two main problems need to be solved: (ⅰ) The liquid-glass transition should be described; we propose such a description in terms of the field theory in the spirit of the fluctuation theory of phase transitions ^[44,45], and we define the field of the "order" parameter corresponding to the glass state of matter; and (ⅱ) the description should be coherent with the PF description of crystallization and imply the competition of these processes.

In our previous work, we proposed for this purpose the topological (or gauge) theory of the glass transition, developed in ^[12,46]. In this theory, the glass transition is a phase transition in a system of topologically protected distortions (disclinations) in some locally ordered atomic structure. This theory has a long history, dating back to the 1970s and 1980s, and it describes the main properties of glass transitions well. The idea of the theory was originally proposed by G. Toulouse ^[47] and has since been extended and refined by other researchers ^{[48,49,50,51,52,53,54,55,56]}. The glass transition in this case can be understood as a collective phenomenon involving the interactions between topological defects in a crystal lattice (the local order should not be necessarily crystalline). These defects are responsible for the structural disorder and mobility of particles in the glassy state, which are the key properties of glasses. This theory provides a useful framework for understanding the glass transition and its relationship to other physical phenomena such as plastic deformation, aging, and mechanical behavior. It describes all the characteristic properties of liquid-glass transitions, and also sheds light on the role of topological defects in other materials, such as polymers and colloids.

2.2.   Model

We consider the liquid within the framework of Jacob Frenkel's conception, which is based on the similarity of the local structure of the liquid (within several coordination spheres) to the structure of the crystal at moderate temperatures ^[57]. This approach explains the similarity of many physical properties of matter in liquid and solid states observed on small spatial and temporal scales ^{[58,59,60,61,62,63,64,65]}. This concept is the basis of the field theories of glass transition, in which the glasses and undercooled liquids are considered as a locally ordered matter with a large number of topologically protected excitations (disclinations) ^{[12,51,55,56]}. In fact, the disclinations act as the main structural elements, and the theory allows one to describe the number of important properties of the glass transition ^[12].

According to the Frenkel's approach, any liquid at the low temperatures has a local atomic order (at the distance of several coordination spheres), similar to that in a crystalline solid, or to an icosahedral noncrystalline ordering. Let one consider the simplest case with a single type of ordering, described by the order parameter ${\bf \Phi }$ ^[66]. When the temperature $T$ is below certain critical value $T < T_c$, the field ${\bf \Phi }$ becomes nonzero in the considered domain, but the long-range order is absent because of the topologically protected excitations (disclinations). Thus, the fluid can be considered as an elastic medium in which both elastic and plastic deformations occur. The presence of the disclinations provides the fluidity of the substance, and elastic deformations determine a free energy $\mathcal{F}$ of the system, which can be written as follows:

where $\lambda$ is the bulk modulus, $\mu$ is the instantaneous shear modulus, and $\hat{\bf u} = u_{ij} = {\mathrm{d} u_j}/{\mathrm{d} x_i} = \nabla_iu_j$ is the distortion tensor (${\bf u}$ is the strain vector). Note that $\mu$ is a microscopic parameter. In the case of the observation times exceeding the short structural relaxation, the averaged (measured) static shear modulus is zero due to the presence of mobile plastic distortions in this system ^[12].

We consider a system in mechanical equilibrium, and ${\bf u}$-field is a free one. The elastic energy is away from the minimum of free energy landscape because the disordered structure is geometrically frustrated. Due to the frustration, ${\bf u}$-field contains stress regions caused by the topologically protected distortions. The topologically protected rotation distortions correspond to the disclinations (or vortex lines) ^[12]. The presence of the disclination breaks the connectivity of the space and leads to the appearance of an irreducible part in the distortion tensor corresponding to the rotation when moving around this disclination:

where the space integration is performed over the dimensionless variable $|{\bf r}| < 1$: $V^{-1}\int\mathrm{d}V = \int\limits\mathrm{d}^3{\bf r}$, and ${\bf \Omega}$ is the Frank pseudovector. The partition function of the system containing $N$ disclinations can be represented in the form of the functional integral:

where $\beta = 1/(k_bT)$, $\delta (\ldots)$ is the functional delta-function, and $J_{\bf r} = \pm 1$ is a topological charge.

The partition function can be written as follows using the functional integral representation of the delta-function:

where ${\bf A}$ is an ancillary field, which follows the condition Eq (2.2), then following ^[12], one can consider the simplification where compressibility was neglected, so the rotation tensor is non-diagonal ($u_{ii} = 0$). After the integration over $\hat{\bf u}$-field, and averaging over the grand canonical ensemble of the vortices (see ^[12], Appendix I), the effective Hamiltonian density assumes the form:

where $g$ is the density of the vortex system. It is nothing else than the Hamiltonian density of the sine-Gordon theory ^[46,67].

Such a topological system as described by Eqs (2.3) and (2.4) undergoes a phase transition ^[12] at the temperature $T_g$:

where $\Omega$ is equal to the Frank pseudovector modulus $\Omega = |{\bf \Omega}|$ ^[12]. Quantum-field theory suggests that, in the three-dimensional case, only the first four terms of the Taylor series expansion of the cosine function with respect to ${\bf A}$ are significant in the fluctuation region near the phase transition point ^[45]. It leads to the following effective Hamiltonian density form:

where $M^2 = g\Omega^2(T-T_g)$ is the square of the effective ${\bf A}$ field "mass".

The model describes a phase transition in a system of interacting topological defects (vortices), leading to the emergence of quasi-long-range order at finite temperature, corresponding to a vortex system with infinite correlation length. In such a model, one considers a topological phase transition, since it is defined by the topological properties of the system and cannot be explained by standard arguments. In ^[66] it was suggested that the $\nabla\times{\bf A}$ vector field could be used as an order parameter for such a transition. In order to define the physical meaning of this field, one follows the liquid state theory of Jakov Frenkel. We assume the presence of a local ordering in the liquid structure. In this case, a long-range translational order is absent if the locally ordered clusters are disoriented with respect to each other. The absence of this disorientation corresponds to $\langle {\bf \hat u }^2\rangle = 0$, which is equivalent to (see ^[66]):

It means the long-range crystalline ordering in the vicinity of the given point. Thus, $|\nabla\times {\bf A}|$ is the degree of distortion of the initial structural order at a given point, and the average of this value indicates the degree of disorder of the structure.

3.   Dynamics

In our previous works ^[66,68,69], we considered some of the simplest cases, where the temperature field was homogeneous, and $T$ was uniform throughout the domain. In the present study, we consider an inhomogeneous temperature distribution in the domain, with a continuous heat flux controlled manually or by solving the heat transfer equation. We assume that the phase transition in the subsystem of topological defects, corresponding to the glass transition, is accompanied by the release of a negligible amount of heat, and the temperature rate change is significantly higher than the rate of the structural changes. Therefore, we neglect the temperature change associated with this transition and focus only on the thermal fluctuations.

3.1.   Equations of motion and heat transfer

Taking into account the presence of thermal fluctuations, the kinetics of the system with a non-conserving order parameter can be described in terms of stochastic nonequilibrium dynamics ^[44,45,70]. The following kinetic equations can be derived from the variational derivative of the full Hamiltonian Eq (2.6) with the addition of the stochastic source as $-\Gamma\partial_t {\bf A} = {\delta \mathcal{H}}/{\delta {\bf A} }+\beta^{-1}\xi$ ^[70]. Thus, the equation describing the dynamics of the ${\bf A}$-field can be written as a stochastic equation taking into account thermal fluctuations in the fluctuation region near the phase transition:

where $\gamma$ is a phenomenological interfacial energy coefficient, and $\xi$ is the dimensionless source of thermal fluctuations, which is the random variable with the normal distribution and unit variance. The local temperature value $T$ obeys the heat transfer equation (or is specified by the given function $T(x, y, t)$ in the degenerate cases described below).

We consider the case of the local thermalization in the representative volumes with size $l\approx50$, which is significantly larger than the size of the nodes $\ell$ exposed to the thermal fluctuations. We believe - that in the case of the absence of feedback from the ${\bf A}$-field equation to $T$, one can introduce the large-scale heat transfer equation as an external field Eq (3.1).

3.2.   Numerical implementation

We have performed numerical simulations of the dynamical Eq (3.1), which was obtained with the Hamiltonian Eq (2.6) in the two-dimensional computational domain with the isolated (zero flux) boundary conditions on ${\bf A}$. In a general case, the three-dimensional case of such a system of equations can also be considered, however, the qualitative results are very similar to those obtained in the two-dimensional study. The symmetry of the system is reduced in Eq (2.3) and the continuous field formulation should not lead to unexpected effects in terms of sensitivity to the cooling kinetics.

The initial distribution was random with $|{\bf A}| = 1$. Such an initial distribution corresponds to the disordered liquid phase. The noise sources were obtained from the random function generator with the normal distribution and zero mean. The computational domains consisted of $L = 300\times 300$ dimensionless units with up to $300$ grid points along the edges; the maximum triangular mesh element size was set as $\ell = 1$. Provided in Table 1, parameters were used in most of the calculations. We tested the mesh and time-step convergence, and for the specific rates there were no significant influence on the size of the topological peculiarities or on the order parameter dynamics during the solidification process. The domain size was also sufficient to proceed to the formation of the bulk multigrain phase.

To compute the given problem, we used the COMSOL Multiphysics 6.0 software ^[71] with PARDISO direct solver and the backward differentiation formula for time integration. The adaptive time step $\delta t = 1.0\times10^{-6}...1.0\times10^{1}$ was controlled at each step on convergence due to the presence of fluctuations and was fixed at the maximum value once the entire domain reached the target temperature $T_1$. All calculations were performed on a dual-processor AMD Epyc 7302-based computer.

3.3.   Quenching simulations

The cold left domain wall simulates the case of ultrafast melt quenching (on the spinning drum), which is the most common experimental setup to produce glass from metallic melts ^[1,7,72,73]. We considered the presented model which experiences bulk quenching with cooling rate equal to $V_{cool} = \partial_t T$ controlled by the heat transfer equation or an arbitrary function $T(x, y, t)$. The initial state corresponds to the disordered liquid with $|{\bf A}| = 1$ at the initial temperature $T = T_{initial} > T_g > T_0$. There are three different cases that have been introduced in the simulations performed:

● Section 4.1: Gradient temperature distribution $T = T_G x$, which was fixed in time during the simulation. No heat transfer equation was implemented in this case.

● Section 4.2: Solution of the heat transfer equation with the Dirichlet boundary condition $T = T_0$ at the left domain wall, $x = 0$, with different thermal conductivity $k$ (see Eq (3.1)). Initial $T$ values were set as $T = T_{initial}$.

● Section 4.3: Initial $T$ was constant $T_{initial}$, then the domain was continuously cooled down with the prescribed $V_{cool} = V_0+V_1x$ at $x = [0..L]$, till $T = T_0$. No heat transfer equation was implemented in this case.

The structural state of the system is characterized by the distortion correlation function $\langle{\bf \hat u\hat u}\rangle\propto \langle (\nabla\times {\bf A})^2 \rangle$ ^[66,69].

To analyze the time correlations in the domain, let one introduce the time-correlation function $G_t[ f({\bf r }, t) ]$ which denotes the simple time correlation function in real space. It describes the spatio-temporal correlations of the given function $f({\bf r}, t)$,

where $\tau_{max}$ is the integration window chosen as half of the computation time $\tau_{max} = t_{max}/2$. In our study, we compute the correlations of the two-dimensional curl of ${\bf A}$-field for $f = (\nabla\times {\bf A})$, which is averaged over all ${\bf r}$ and corresponds to ${\bf k} = 0$ in the reciprocal space.

4.   Results and discussion

4.1.   Fixed gradient temperature distribution

The simulations performed (see Figure 1) allow one to define the structure selection depending on the final temperature distribution. The figures represent the gradual coarsening of the structure, which is different depending on the position in the domain and the corresponding $T$. Qualitatively similar structure selections depending on the different final temperatures performed at different cooling rates $V_{cool}$ can be found in ^[74].

For the temperatures $T < 0.02$, the amorphous phase appears on the left side of the domain $x \in [0..20]$ in the cold region, containing the structural inhomogeneities with the size around $q \approx 0.2$. This region is preserved at the very cold temperatures and correspondingly low fluctuations, which cannot have any observed effect on the further coarsening of the structure. Following to the right side, to the higher temperature region, the transition zone can be found at $x \in [20..40]$, which is characterized by the varying inhomogeneities of size $q \approx 5$ corresponding to the relatively stable fine-grained amorphous phase, which stays at $T < 0.04$ and experiences some coarsening up to $q \approx 8$, conditioned mainly by the ratio of surface and volume energy, resulting in the relatively stable microstructure. The fine polycrystalline structure can be observed from $x \in [40..120]$ at the later stages of the structure relaxation. This structure is formed from the coarse amorphous phase (can be seen in Figure 1(a)) during the annealing, which is mainly driven by the stochastic fluctuations of the ${\bf A}$-field. The central zone $x \in [120..180]$ is filled with the crystalline phase with the relatively large size of grains $q\approx 50$ growing from the coarse polycrystalline structure (see Figure 1(a)) with multiple nuclei. These oriented crystals contain a tiny amount of defects. The formation of the large crystals occurs mostly at $t \in 2500..5000$ after the primary coarsening. The phase boundary is located at $x = 187$ at $T = T_g = 0.187$, which is the equilibrium transition temperature $T_g$. The phase to the right of the $x > 187$ region corresponds to the liquid at high $T > T_g$.

4.2.   Influence of the finite thermal conductivity coefficient

To study the influence of the cooling with the varying heat flux controlled by the unsteady heat transfer equation, Eq (3.1), on the final microstructure, we performed a series of simulations with different heat transfer coefficients $k$. The Dirichlet boundary condition $T = T_{0}$ applied to the left domain wall allows one to model the case where the cold wall is in contact with the cooling reservoir as in rapid melt spinning experiments with the large cooled spinning drum. The thermal conductivity $k$ here controls the total heat flow through the melt.

4.2.1.   Effect of shear modulus

One can distinguish between the formation of the columnar crystals with the cross-section size $Q\approx 45$ in the slow cooling regime, at $k = 5$, which is presented in Figure 2. The left domain wall has cooled, becoming a region of the heterogeneous nucleation. The crystal boundary follows the isotherm $T = 0.165$, which corresponds to the effective temperature of the phase transition in the $A$-subsystem. The resulting temperature on the front is lower than the equilibrium one $T_g = 0.187$ observed in the case of the fixed temperature gradient (see Section 4.1, Figure 1). During the later cooling (see Figure 2(b)), the structure shows the result of the homogeneous nucleation and sequential polycrystalline growth. The ordered islands of the significant size $q\approx 32-45$ are comparable to the cross-section of columnar crystallites formed during the slow cooling and can be considered as equiaxed crystals forming polycrystalline aggregates. In our case, such polycrystalline structure resulted from the bulk relaxation on the later stages of the cooling. During the process of microstructure formation, the following stages are observed: (ⅰ) Emergence of fine amorphous structure after cooling below $T_g$, (ⅱ) coarsening of the fine structure under the conditions of relatively high thermal fluctuations (corresponding to the annealing at $\approx 0.7-0.8 T_g$), and (ⅲ) ordering and consolidation of adjacent regions and migration of defects to the grain boundaries.

The $T_g$ isotherm accelerates to the right side of the domain, leading to the columnar to equiaxed transition (similar to that observed during the crystallization in the many metallic alloys ^[75,76]). The effect of the limitation of the crystallization front propagation velocity is caused by the finite velocity of the reconfiguration of the ${\bf A}$-field, which is controlled by the $\Gamma_{\bf A}$ coefficient. In the simulations here we found that the transition to the bulk crystallization occurs at $t\approx8000$ and the front position up to this time can be interpolated by the 2$^{nd}$-order polynomial: $X = 7\times10^{-7}t^2+2.7\times10^{-3}t+35.375$, so the maximum front propagation velocity is $V_{max} = 0.0139$ dimensionless units.

Note that $\mu$ is a microscopic parameter that is not zero in the liquid state. It corresponds to the macroscopic shear modulus in the case when the relaxation time of the structural system significantly exceeds the observation time. Usually, the shear modulus is inversely related to the temperature. Here, in our model $\mu$ controls the strain (inverse surface energy) of the ${\bf A}$-field, and for the convenience we assign the additional parameter $\gamma$ for the precise selection of the ratio of the equilibrium transition temperature to the surface energy generated by the inhomogeneities of ${\bf A}$. To demonstrate the effect of the increasing $\gamma = 1$, one can refer to Figure 3, where for the low diffusivity coefficients $k = 1, 5$, the formation of the elongated crystals stops at once. The shear modulus leads to the disordered crystal growth Figure 3(a). Comparing the columnar growth with the same $k = 5$ from Figure 2(a) to Figure 3(b), there's a huge difference: Although the $T_g$ and the kinetic undercooling are the same, the structure is significantly grained. The shear modulus value directly influences the minimum possible grinding curvature as a surface energy factor.

4.2.2.   Effect of the thermal conductivity coefficient

Figure 4 shows the intermediate stages and the final structures after the cooling with $k = 10$ (Figure 4(a), (b)), and $k = 20$ (Figure 4(c), (d)). Due to the simple $T$ distribution, which is very close to the $erfc()$ solution of the heat transfer equation ^[77] in semi-infinite media, the right side of the domain experiences slower cooling rates in the later stages. However, the intermediate zone of the domain, $x \in [30..150]$, is subject to longer heat exposure and, hence, structural relaxation. One can distinguish (see Figure 4(a)) between the formation of the large polycrystalline structure with the sizes of around $q \approx 32$. Here, the local cooling rates are comparable to those in the $k = 5$ case. However, the right panel Figure 4(b) shows the slight coarsening of the structure in the intermediate region $x \in [30..150]$ up to $q \approx 35$. The gradient of the properties in the case $k = 10$, Figure 4(b) is caused by the long heat exposure (when $T\in [T_0..T_g]$) in the intermediate zone before the final freezing of the structure. The left panel of Figure 4(c), $k = 20$, shows the propagation of fine amorphous structure $q\approx 7$ (transition zone is blurred) into the liquid, which is relatively fast cooled, and thus forms the relaxed amorphous structure of $q\approx 20$ shown in Figure 4(d). One can still find the narrow zone $x\in[50..130]$ of the coarsened intermediate structure in Figure 4(d) with $q \approx 28$, but its boundaries are blurred. The effect of the crystallization front drag and transition to the vitrification (similar to that observed during the solidification of Cu-Zr ^[34,35]) here can be observed here on the small grains on the left side of the domain Figure 4(c), which later appeared coarsened in the Figure 4(d) due to the subsequent slow relaxation. The interesting feature of the final structures in Figure 4(b), (d) is the presence of the distinguished elongated defect chains.

The almost instantaneous solidification (finished in the whole domain in $\Delta t\approx 1750$) occurred during the cooling with the thermal conductivity $k = 50$, see Figure 5. In Figure 5(a), the front of the amorphous phase can be clearly distinguished, while the phase front is rather blurred in present case ($\delta Z \approx 45$). After the transition to the $|{\bf A}| = const$ phase a gradual relaxation with coarsening from $q = 17$ to $q = 23$ was observed. After the relative long relaxation (around $\Delta t = 25000$), no further relaxation was observed. This can also be shown with the time correlation function $G_t(\nabla\times {\bf A})$ (see Figure 6). It should be noted that the final amorphous state of the disclination subsystem contains clusters that coalesce into a low-dimensional foam-like structure. The formation of such structures has often been observed in computer models of glasses (see, for example, ^[78]).

Obtained distributions of structures from Figures 2–5 qualitatively correspond to the drums spinning at different speeds, and, thus, cooling rates ^[79] for Fe-Si-B ^[80], rare-element based alloys ^[81], Al-Si ^[82]. The columnar to equiaxed transition, Figure 2, is observed in spinning drum experiments with different materials such as Cu-P-Ni ^[83] and Co-based alloys ^[84], and the gradient properties were observed in experiments with Ni-based alloys ^[85] and coarsening in ^[86]. The effect of cooling rates are well studied for Fe-based ribbons ^[80,87] and rare-element based alloys ^[81], and are also qualitatively reproduced by our model.

4.2.3.   Time correlation function in subdomains

The time correlation functions $G_t(\nabla\times {\bf A})$, Eq (3.2), calculated for $t\in[0..10000]$, show different slopes of $G_t$ for the amorphous and crystalline phases in the subdomains; see Figure 6. Here, the upper curve corresponds to the amorphous phase formed for $k = 50$; see Figure 5(c). The amorphous curve represents two different stages of relaxation (short-time and long-time) and crossover between them with shallow slope of the long-time correlations. The transition curve for $k = 10$ (see Figure 4(b)) represents the amorphous-like correlation in the fine-grained polycrystalline zone. The left part of the slowly cooled domain $k = 5$ with the oriented crystals shows the three-stage relaxation and a decrease of correlations on long timescales (see Figure 2). The right side of the domain $k = 5$ contains the liquid disordered phase and, hence, the absence of correlations. The total short-time correlations show an exponential relaxation associated to the mobility of the ${\bf A}$-field during primary cooling at $t\in[0..1000]$ in all solidified domains. The relaxation in the amorphized subdomain, Figure 6, is flatter due to the slowing down of the dynamics since the domain is at low $T$ most of the time. The liquid disordered phase shows almost instantaneous relaxation.

4.3.   Controlled cooling speed

The simulation with finite thermal conductivity values introduced earlier in Section 4.2 shows that the cooling rate of a local sample region depends on time and its spatial coordinates. In the current section, we have selected the most representative cooling rates and set the cooling rate at each point of the domain as $V_{cool}(x) = V_0+V_1x$ at $x = [0..L]$. The initial temperature was set as $T = T_{initial}$. The resulting microstructure of the solidified domain is shown in Figure 7(a). Here, one can see a pronounced effect of the $V_{cool}$ on the final relaxed microstructure of the solidified domain. The cooling rate has a significant effect on the structure formation. Here, we obtain the gradient of simultaneously selected structures: From columnar crystals to polycrystalline and then to the amorphous structure, see Figure 7(a). With the cooling rate gradient, the whole range of the formed structures can coexist with each other at the low temperature (since they correspond to a single minimum of the free energy of the gauge parameter ${\bf A}$). The carried out numerical analysis of the kinetics of our model under the conditions of rapid quenching shows the qualitative agreement with the microstructures formed during rapid quenching ^{[86,88,89,90]}.

A similar sensitivity to the cooling rate has been observed in previous studies of the gauge models of vitrification ^[66,68,69], although much smaller domains with full thermalization were considered. The distinction between crystalline and amorphous phases as final structures lies in the path to the final structure, which can be summarized as the evolution of the order parameter $\langle (\nabla \times {\bf A})^2 \rangle$ ^[66]. Previously, a hypothesis was proposed that the presence of a maximum on the $\langle (\nabla \times {\bf A})^2 \rangle_t$ evolution curve distinguished the presence of the amorphous state ^[66]. We assume that the glass transition occurs when the first derivative of the order parameter becomes zero as a result of the drastic change in kinetics from the critical to the relaxation regime. After the sharp increase up to the maximum value resulting from the glass transition, the order parameter for glass starts to decrease slowly as a result of aging. These differences are clearly visible in Figure 7(b), where x = [0..50] corresponds to the crystallized subdomain, x = [50..100] to the mixed/polycrystalline zone, and all other domains are vitrified. The higher value of the glass order parameter corresponds to a more amorphous system. The low magnitude of the averaged $\langle(\nabla\times{\bf A})^2\rangle_t$ on the plateau (for later times, $t > 5000$) represents the presence of a small number of inhomogeneities and defects in the volume, which is typical for crystals with isolated defects or individual grains with pronounced boundaries.

5.   Conclusions

The proposed model describes the solidification of the melt during the quenching simulations with a cold wall boundary condition, which is assumed to be in contact with the cooling reservoir as in rapid melt spinning experiments with the cooled drum. In the early stage, the solidification process starts with the formation of a fine-grained short-range ordering. Such local ordering is fragmented by the topologically protected excitations (disclinations), which prevent the formation of long-range translational order. The phase transition in the system of topologically protected excitations, which is implemented as the vector gauge field ${\bf A}$, was found to be sensitive to the local cooling rates, which are directly related to heat fluxes and heat thermal conductivity coefficients.

Numerical simulations yielded the cooling rates corresponding to the competing regimes of crystal growth or amorphous phase formation. In the case of a high thermal conductivity $k = > 20$ and, hence, high cooling rates, the domain solidifies into a bulk amorphous state with the emerging blurred transition zone revealing high magnitudes of the corresponding order parameter $\langle (\nabla \times {\bf A})^2 \rangle$. The structure corresponding to this state consists of a system of topological defects that tend to form low-dimensional clusters, forming a complex and disordered network. As the thermal conductivity of the material decreases $k < = 5$, a boundary between the polycrystalline and columnar crystal appears. The largest ordered (crystalline) columnar aggregate, where the field $(\nabla \times {\bf A})^2$ magnitude is small, has been observed for the case of slow solid-liquid front propagation. With the increased cooling rate and the corresponding increase of driving forces in the transition zone, we observed a columnar to equiaxed transition (CET) qualitatively similar to that observed in many metallic alloys ^{[75,76,83,84]}. The moving front follows the equilibrium isotherm corresponding to the transition temperature $T_g$, although front drag leads to the effect of kinetic undercooling and the emergence of the crystallization front maximum velocity.

With the cooling rate gradient, a whole range of the formed structures can coexist with each other at the low temperature (since they correspond to a single minimum of the free energy of the gauge parameter ${\bf A}$). Nevertheless, even the fine-grained amorphous structure gradually coarsened up to the certain value $q \approx 23$, which was revealed during long simulations and corresponds to the magnitude of the parameter $\mu^{-1}\beta^{-2}$.

Despite the fact that we considered the simplest case, with a single ordering type, and corresponding topologically stable excitations, the model exhibited main characteristics typical for a glass-forming alloys. In particular, the behavior of this model showed competition between the vitrification and crystallization processes. The carried out numerical analysis of the kinetics of our model under the conditions of rapid quenching with simultaneous solution of the heat transfer and ${\bf A}$-field equations shows the qualitative agreement with the microstructures formed during rapid quenching in rapid melt spinning experiments ^{[79,80,81,82,85,86,87]}.

Appendix I

The effective Hamiltonian of our system is written as follows:

Carrying out the averaging over a grand canonical ensemble of the "particles" endowed with the two possible dimensionless charges $J_n = \pm 1$, there are possible vortices' configurations that can be taken into account. Then the path integral is as follows:

After summation, the effective Hamiltonian density takes the following form:

where $g = e^{-\beta\mathcal{E}_c}$ is the density of vortices.

If one expands the cosine term into the power series, taking into account that in the 3D case only the first two terms of this expansion are relevant ^[12], then one obtains the following effective form of the Hamiltonian density:

One can represent the ${\bf A}$ field as the sum of fast, $\tilde{\bf A}$, and the slow, ${\bf A}$, parts if close to the critical point: ${\bf A}\to {\bf A}+\tilde{\bf A}$. After averaging over the fast field, the Hamiltonian density is written as follows:

where

The shear modulus usually inversly depends on $T$. We assume that near the glass transition this dependence is a linear function: $\mu = \varepsilon (T_{\mu }-T)$, where $T_{\mu }$ is some temperature at which the shear elasticity appears in the liquid. Taking this into account, we lead to the following effective Hamiltonian density presentation:

where $M^2 = g\Omega^4\varepsilon (T-T_g)/(2\pi)^2$ is the square of the effective ${\bf A}$-field "mass", and

is the phase transition temperature in the disclination subsystem. From above, the glass transition temperature is proportional to the shear modulus of the glass state, $T_g = \mu(T_g) (\Omega/2\pi)^2/k_b$. Thus, it is the well-determined thermodynamic value directly related to the elastic properties.

Author contributions

V. Ankudinov: Conceptualization, Data curation, Formal analysis, Investigation, Visualization, Writing-original draft, Writing-review & editing, Formal analysis, Validation; M. Vasin: Conceptualization, Data curation, Investigation, Writing-review & editing, Formal analysis, Funding acquisition, Supervision; K. Shklyaev: Investigation, Validation, Writing-review & editing, Visualization. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The work was supported by the Russian Science Foundation under the project No. 23-22-00168.

Conflict of interest

No potential conflict of interest was reported by the authors.