Dynamics of grain boundary premelting

The mechanical strength of a polycrystalline material can be drastically weakened by a phenomenon known as grain boundary (GB) premelting that takes place, owing to the so-called disjoining potential, when the dry GB free energy \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{gb}$$\end{document}σgb exceeds twice the free energy of the solid–liquid interface \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{sl}$$\end{document}σsl. While previous studies of GB premelting are all limited to equilibrium conditions, we use a multi-phase field model to analyze premelting dynamics by simulating the steady-state growth of a liquid layer along a dry GB in an insulated channel and the evolution of a pre-melted polycrystalline microstructure. In both cases, our results reveal the crucial influence of the disjoining potential. A dry GB transforms into a pre-melted state for a grain-size-dependent temperature interval around \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_m$$\end{document}Tm, such that a critical overheating of the dry GBs over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_m$$\end{document}Tm should be exceeded for the classical melting process to take place, the liquid layer to achieve a macroscopic width, and the disjoining potential to vanish. Our simulations suggest a steady-state velocity for this transformation proportional to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{gb} -2 \sigma _{sl}$$\end{document}σgb-2σsl. Concerning the poly-crystalline evolution, we find unusual grain morphologies and dynamics, deriving from the existence of a pre-melted polycrystalline equilibrium that we evidence. We are then able to identify the regime in which, due to the separation of the involved length scales, the dynamics corresponds to the same curvature-driven dynamics as for dry GBs, but with enhanced mobility.


Dynamics of grain boundary premelting M. Torabi Rad, G. Boussinot * & M. Apel
The mechanical strength of a polycrystalline material can be drastically weakened by a phenomenon known as grain boundary (GB) premelting that takes place, owing to the so-called disjoining potential, when the dry GB free energy σ gb exceeds twice the free energy of the solid-liquid interface σ sl . While previous studies of GB premelting are all limited to equilibrium conditions, we use a multi-phase field model to analyze premelting dynamics by simulating the steady-state growth of a liquid layer along a dry GB in an insulated channel and the evolution of a pre-melted polycrystalline microstructure. In both cases, our results reveal the crucial influence of the disjoining potential. A dry GB transforms into a pre-melted state for a grain-size-dependent temperature interval around T m , such that a critical overheating of the dry GBs over T m should be exceeded for the classical melting process to take place, the liquid layer to achieve a macroscopic width, and the disjoining potential to vanish. Our simulations suggest a steady-state velocity for this transformation proportional to σ gb − 2σ sl . Concerning the polycrystalline evolution, we find unusual grain morphologies and dynamics, deriving from the existence of a pre-melted polycrystalline equilibrium that we evidence. We are then able to identify the regime in which, due to the separation of the involved length scales, the dynamics corresponds to the same curvature-driven dynamics as for dry GBs, but with enhanced mobility.
In a polycrystalline material at a temperature well below the bulk melting temperature T m , the disordered region at a Grain Boundary (GB) is only a few atomic distances wide. As the temperature increases towards T m , however, GBs may turn into a liquid-like, thermodynamically stable thin (i.e., nanoscale) film. This order-disorder transition is termed GB premelting 1 . Premelting may be observed not only at the GBs but also on the surfaces 2,3 and, due to its physically-rich nature and important consequences in wide range of applications, has interested scientists from not only chemistry and physics domains, but also from domains such as earth science 4 , fluid mechanics 5 , and, interestingly, biology 6 . GB premelting has implications in, for example, general metallic systems 1 , steel 7 , and ceramics 8 . It increases GB diffusivity and mobility 9,10 , and can drastically influence the macroscopic properties of a material. For example, it can reduce the resistance to shear stresses 11,12 , which can result in material failure in high-temperature applications. GB premelting is also linked to cracking during late-stage solidification of metallic alloys 13,14 , to liquid metal embrittlement 15 and in general influences grain growth, for example in situations such as brazing. Despite these implications, GB premelting is still not fully understood one reason being the lack of abundant experimental studies (see 16 and references therein), which is mainly due to the challenges in observing internal material interfaces such as GBs. Another reason is that the existing experimental evidence for pure materials is controversial 12 due to a fair criticism that the conclusions could have been influenced by trace impurities 17 .
In the literature, basic thermodynamics are invoked to suggest that pre-melted GBs form because solid-liquid interfaces may repulse each other and form a liquid-like layer to decrease the free energy of the system 12 . Generically, the energy σ gb of a low-angle GB is small and proportional to the misorientation between the adjacent grains 18 , while σ gb is higher for high-angle GBs. Moreover, σ gb may also be influenced by the presence of GB adsorption. Then, when σ gb is larger than twice the energy of well-separated solid-liquid interfaces σ sl , the system can lower its free energy by replacing the dry GB with two solid-liquid interfaces and a liquid layer between them. This is quite similar to the situation where a solid in contact with vapor starts to melt at its surface when it is energetically favorable to replace the solid/vapor interface by a solid/liquid interface and a liquid/ vapor interface 19 .
The thickness of the liquid layer at a pre-melted GB is determined by the competition, as the distance between the solid-liquid interfaces varies, between the changes in bulk and interfacial free energies. Theoretically, the variation of interface energy is provided by a so-called disjoining potential, that describes energetically the deviation of the atomic structure in the thin liquid layer from the one in the bulk liquid and the interaction between the solid-liquid interfaces. This potential allows the interface energy to be interpolated continuously between σ gb , at vanishing distance between the solid-liquid interfaces, and 2σ sl at infinite distance 20 .
Several computational methods were employed in order to study premelting. Using Lennard-Jones or embedded-atoms-method potentials, Molecular Dynamics [21][22][23] and Monte-Carlo simulations 24,25 were performed on the atomistic level. Phase-field or diffuse-interface methods were also used. In Refs. [26][27][28] , one order parameter discriminates between solid and liquid, and one orientation field describes crystallographic orientation in the solid phase. In Refs. 16,29,30 , a multi-phase-field model is used, in which the liquid and each grain are characterized by a separate field. More recently, the phase-field-crystal method was also employed 12,[31][32][33] , giving also some hints on the microscopic mechanisms of premelting.
In general, premelting under equilibrium conditions is mainly controlled by the dependence of the disjoining potential on the distance between the two solid-liquid interfaces and by the variation of this dependence with temperature. In the simplest case, the disjoining potential is monotonous for any temperature, and the equilibrium width of the liquid layer varies continuously with the temperature. As soon as the disjoining potential becomes non-monotonous, the scenario becomes more complicated. In particular, so-called "thin-to-thick" first-order transitions may occur when the disjoining potential corresponds to damped oscillations 25 , with a coexistence in a certain temperature interval of several equilibrium widths.
The above studies have resulted in valuable insights into the physically-rich nature of GB premelting. They have, however, to the best of our knowledge, considered only premelting under equilibrium conditions. In this article, we are interested in the dynamics of premelting, i.e., how a liquid layer penetrates along dry GBs driving the transformation of the latter into a pre-melted equilibrium. More generally, we consider the situation where a fully-solid and low-temperature polycrystalline structure is brought to a higher temperature close to T m . Fundamental open questions include the effect of the grain size on the dynamical behavior of a pre-melted layer, the dependence of the penetration velocity on the GB and solid/liquid interface energies, the shape of the liquid layer around its tip, and the influence of the initial temperature on the final structure, i.e. the product of the transformation. We address those questions using the multi-phase-field model with obstacle potentials that is implemented in MICRESS 34 . While the phase-field-crystal model and molecular dynamics techniques have been used for studying premelting in equilibrium conditions, the phase field model is used here to resolve the length scale associated with the diffusion field, too large to be addressed with the former methods. We first present the solution for the pre-melted equilibrium that this model yields. Assuming a steady-state growth of the liquid layer along the dry GB, we analyze the consequences of energy conservation on the equilibrium state that develops as a product of the transformation. We then present simulations that focus on the triple junction and on the dependence of the growth velocity on the interface energies. Next, we present the simulation of the evolution of a polycrystalline structure.

Theoretical phase-field analysis
Phase field models describe the evolution of continuous fields that may represent for example the local temperature, chemical composition or fraction of each phase under consideration. In the latter case, the fields are called 'phase fields' and their variations represent interfaces, for example between two grains or between a solid and a liquid. They may also represent triple junctions when three fields are involved, or even higher order junctions.
The premelting of a GB is linked to the existence of a so-called disjoining potential between solid/liquid interfaces that are separated by small enough distances. This interaction between solid/liquid interfaces describes the fact that the surface energy may be considered as going continuously from the value of the GB energy σ gb at vanishing distance (i.e. when the two solids grains are in contact), to twice the value of the solid/liquid interface energy σ sl at infinite distance. A certain variety of premelting scenarios can be expected depending on the shape of the disjoining potential 27,28 . When the latter varies monotonically with the distance, the scenario is however rather simple and premelting occurs when σ gb − 2σ sl > 0 , with a unique solution linking the width of the liquid layer to the temperature at equilibrium.
In the phase field model, the disjoining potential is related to the overlap of the phase fields of the two grains. There are mainly two ways of describing the dynamics of the phase fields. In the first one with so-called double well potentials, the phase fields exhibits spatial variations of an hyperbolic tangent type, with an exponential convergence to the value 0 or 1. In this frame, it was shown 29 that premelting is not properly reproduced, especially because, although being able to be repulsive at short distances, the disjoining potential is always attractive at large ones, regardless of the ratio of GB to solid/liquid energies. In the second one, which we use in this article, with so-called obstacle potentials, the phase fields spatial variations are strictly restricted to a finite length scale and the fields reach exactly their value 0 or 1 at a finite distance from the center of the interface. The corresponding evolution equations are presented in the "Methods". It was shown (see Fig. 10 in Ref. 30 ) that such a model provides, regardless of the distance, a repulsive (attractive) interaction when σ gb − 2σ sl > 0 (< 0) . This model thus reproduces premelting in a manner quite similar to a sharp interface description using a monotonous disjoining potential, for example a decaying exponential as in Ref. 35 .
Premelting equilibrium. The pre-melted equilibrium reproduced by the multi-phase field model with obstacle potentials, was analyzed in Ref. 30 , and we give here a brief overview of that analysis. The problem is one-dimensional and in the following all lengths are expressed in units of η/(2π) , where η is the width of the interface, an intrinsic property of the phase field model. Three phase-fields φ 1 (x), φ 2 (x) and φ 3 (x) are defined so as to describe the existence of three phases: two solid grains (solid 1 with φ 1 and solid 2 with φ 2 ) and the liquid ( φ 3 ). The sum of these phase fields equals unity at any position in space: The problem is invariant in exchanging solid 1 and solid 2, so that the fields obey some symmetry when x = 0 is chosen as the center of the liquid film: www.nature.com/scientificreports/ Three regions may then be defined, depending on the local number of phase fields that assume a non-vanishing value, as depicted in Fig. 1 in different colors. The black solid line represents φ 1 , the black dashed line represents φ 2 and the blue line represents φ 3 . Corresponding to region (I), we have φ 1 (−∞ < x < −x 1 ) = 1 and φ 2 (x 1 < x < +∞) = 1 . Corresponding to region (II), we have φ 3 (−x 1 < x < −x 2 ) = 1 − φ 1 (x) , which means that only solid 1 and the liquid are present at −x 1 < x < −x 2 ( φ 2 = 0 ), and only solid 2 and the liquid are present at x 2 < x < x 1 ( φ 1 = 0 ). In region (III), all three phase fields are non-vanishing. It is important to note that It can be seen within this frame that φ 3 may not reach unity. We are thus in the case of a pre-melted system for which a liquid layer presenting a certain degree of crystalline order ( φ 3 = 1 ) separates two solid grains. The larger is φ 3 at the center of the liquid film, i.e. φ 3 (x = 0) , the larger is the disorder of this pre-melted GB. When x 2 increases (and x 1 decreases consequently owing to x 1 + x 2 = 2π ), the level of disorder in the liquid layer decreases ( φ 3 decreases), and when x 2 = x 1 = π , φ 3 vanishes identically throughout the whole system, i.e. the GB is dry. On the other hand, when x 2 decreases, i.e. x 1 increases, the liquid phase field φ 3 increases, and when x 2 = 0 , i.e. x 1 = 2π , the level of disorder at the center of the liquid layer is characteristic of the bulk liquid, i.e. φ 3 (x = 0) = 1 . When x 2 decreases further, i.e. when x 2 < 0 , region (III) disappears, the two solids are no more interacting via the disjoining potential ( φ 1 and φ 2 are no more overlapping), and a macroscopic liquid film having the structure of the bulk liquid, with φ 3 (x 2 < x < −x 2 ) = 1 , is then in equilibrium with the two solids. This situation is depicted in Fig. 2.
In the following we denote as a pre-melted GB the situation where 0 < φ 3 (0) < 1 and a melted GB when φ 3 (0) = 1 . Thus premelting occurs in a finite range of x 2 , i.e. 0 < x 2 < π . This is in contrast to a sharp interface description where the disjoining potential exhibits a exponential decay with the distance between the solid-liquid interfaces. Indeed, in this case, the magnitude of the disjoining forces never strictly vanishes, and the width of the liquid layer at equilibrium diverges logarithmically when the temperature T approaches from below the bulk melting temperature T m . Here, in the phase field model, 0 < T m − T → 0 means 0 < x 2 → 0 , and x 2 < 0 corresponds to T = T m .
The analytical solution of the phase field equations for the pre-melted equilibrium, i.e. for 0 < x 2 < π , reads as follows (see Ref. 30 for more details): Phase fields profiles for the pre-melted equilibrium corresponding to a liquid layer, represented by φ 3 (x) , between grain 1, represented by φ 1 (x) , and grain 2, represented by φ 2 (x) . The sum rule i φ i (x) = 1 holds. In region (I) at negative (positive) x, grain 1 (grain 2) is in its thermodynamically stable bulk state with φ 1 = 1 ( φ 2 = 1 ). The liquid layer is present with φ 3 (x) = 0 in regions (II) and (III) (see text for more details). The degree of disorder in the liquid layer φ 3 is characteristic of premelting with φ 3 < 1. www.nature.com/scientificreports/ and with and L corresponds to the latent heat of fusion related to the entropy difference between solid and liquid at T m . The dimensionless quantity σ * = σ gb /σ sl is a measure of the tendency for the GB to spontaneously melt and premelting occurs when σ * > 2 . It can be seen that within this model σ * should remain smaller than 4 for premelting to be described in the terms presented above. As mentioned above, the level of disorder of the pre-melted GB may be quantified by For a given , while x 1 = 2π − x 2 (owing to x 1 + x 2 = 2π ), x 2 is found as the solution of This explicit relation between the temperature and x 2 , that was not given in Ref. 30 , will be used in the following paragraph where we study energy conservation during steady-state propagation of the liquid film along the GB.
Let us now consider the two transitions that are described above. For x 2 = x 1 = π , corresponding to φ 3 (x) = 0 , we have � = −(σ * − 2)/2 . This temperature corresponds to the onset of premelting, i.e. it is the minimum temperature at which premelting can be expected. Note that in the neighborhood of this transition, i.e. www.nature.com/scientificreports/ when x 1 − x 2 ≪ 1 , numerics are tedious since the grid spacing should much smaller than x 1 − x 2 . For x 2 = 0 , one has = 0 , corresponding to the melting temperature, i.e. the onset of melting. The level of disorder at the center of the liquid film is then characteristic of the bulk liquid, i.e. φ 3 (0) = 1 . It should be noted that the onset of premelting approaches the onset of melting when σ * approaches 2. Then, x 2 changes from π to 0 and φ 3 (0) changes from 0 to 1, when changes within a small interval from −(σ * − 2)/2 to 0. When x 2 < 0 as depicted in Fig. 2, the solution for φ 1 in the range −2π + x 2 < x < x 2 (i.e. for the interface between grain 1 and the liquid) is This is the solution for an interface in a two-phases system as for example during dendritic solidification.
Dynamics of GB premelting. As mentioned in the Introduction, the main objective of this paper is studying the premelting dynamics. We consider a set up presented in Fig. 3 where the liquid layer propagates at a constant velocity V along the dry grain 1/grain 2 boundary in a two-dimensional system with a width . The growth is accompanied by heat diffusion and the system is insulated. Thus no heat flux is assumed at the lower and upper boundaries of the box, i.e. in the direction normal to the growth. The temperature T ∞ of the dry GB is represented by the dimensionless quantity which is a given of the problem, and which is the temperature that is prescribed at the right boundary of the box. During steady-state growth, as we will see later, heat fluxes are present in a restricted region around the tip of the liquid layer, and, far enough behind, the system equilibrates at a temperature T −∞ represented by and, in general, −∞ � = ∞ . The equilibration far behind the tip leads to the one-dimensional equilibrium described in the previous section, with −∞ playing the role of .
This set-up is typical for the study of a confined growth regime, for example the dendritic growth in a channel. The new phase, here the liquid, is assumed to nucleate at the GB and subsequently the liquid layer enters the growth regime. This scenario is characteristic of a first-order phase transition for which thermal fluctuations are needed and an energy barrier has to be overcome in order to locally produce a seed of the new phase. While premelting at free surfaces does not seem to require such an energy barrier and takes place homogeneously over the surface 19 , it should be different in the case of a grain boundary since an elastic accommodation is in general required due to the difference in atomic density between the solid and liquid states 36 . We thus assume that nucleation events are rare enough and the nucleii are far enough from each other, so that the liquid layer is able to grow steadily along the GB. Therefore, our set-up corresponds to a situation where the polycrystalline structure, stable at low temperature, is suddenly brought to a higher temperature ∞ . In an experimental situation, a time is needed for the system to reach this temperature; hence, our set-up is relevant to experiments when the time for nucleation is larger than the time for thermal equilibration at ∞ . Note that such a thermal equilibration may take place on a very short time scale in a thin-sample with two-dimensional dynamics as it is the case here.
The description of the pre-melted equilibrium in the previous section allows to find, for a given temperature, the geometrical characteristics of the liquid layer, i.e. x 1 and x 2 . In the case of a growth of the liquid layer at constant velocity V, conservation of energy implies that the energy far ahead of the tip equals the energy averaged over the direction perpendicular to growth far behind the tip where equilibration takes place. In other words, the enthalpy associated with the creation of the liquid is compensated by changes in temperature and interface energy. This yields an equation that relates −∞ to ∞ , through the equilibrium profiles φ i (x) given in the previous subsection ( i = 1, 2, 3 ) and their spatial derivatives φ ′ i (x): Figure 3. Schematics of the two-dimensional steady-state growth at velocity V of the liquid layer along the dry boundary between grain 1 and grain 2. Ahead of the liquid's tip, the prescribed dimensionless temperature is ∞ , and far behind the tip, the system equilibrates at a dimensionless temperature −∞ .

Scientific Reports
| (2020) 10:21074 | https://doi.org/10.1038/s41598-020-77863-9 www.nature.com/scientificreports/ where L = ηL/σ sl and S = 4c p T m /L with c p the specific heat. The left hand-side of Eq. (15) corresponds to the energy far ahead of the tip, while the right hand-side corresponds to the energy far behind the tip. This equation may be solved numerically by finding which −∞ , associated with its corresponding x 1 and x 2 , is fulfilling the condition (15). It can be seen that the solution for −∞ depends on the degree of confinement . In Figs. 4 and 5, we present the dependence on of −∞ and x 2 respectively, for three different values of σ * > 2 , and for ∞ = −0.047 , S = 40 and L = 15 . These values for S and L , that can correspond for example to T m = 800 K, and a characteristic width of the liquid layer η = 7.5 nm, are typical for metals. In particular, they are close to the properties of pure Al, for which some experimental results are compared to our polycrystalline simulation later on. In Figs. 4 and 5, it can be seen that −∞ increases with and x 2 decreases when increases, corresponding to an increases in the width of the liquid layer when increases. At large , −∞ converges for all σ * to the dashed line that corresponds to ∞ . The limit � −∞ ( → ∞) → � ∞ results from the fact that the terms proportional to then dominate in Eq. (15), the energy increase caused by the presence of the liquid layer (responsible for −∞ ≤ ∞ ) becoming negligible in the energy conservation relation. Correspondingly, the geometry of the liquid layer represented by x 2 in Fig. 5 becomes independent of when ≫ 1 , and may be estimated replacing by ∞ in Eqs. (10) and (11).
On the other hand, when becomes small, it can be seen from Fig. 5 that x 2 goes to π . As mentioned in the previous paragraph, this situation corresponds to the onset of premelting, and indeed � −∞ ( → 0) → −(σ * − 2)/2 in Fig. 4. Let us note however that may not be arbitrarily small since it should, at least, exceed the interface width, i.e. 2 π in our dimensionless units.
The -dependence of −∞ and x 2 implies a -dependence of the level of disorder of the liquid layer that grows along the dry GB. As can be indeed seen in Fig. 6, while the level of disorder φ 3 (0) given in Eq. (9) intuitively increases with σ * , it also decreases when decreases (since −∞ increases with ).
Let us note that, as a consequence of the fact that ∞ is close to −(σ * − 2)/2 for σ * = 2.1 , the variations of −∞ and x 2 with are small. Thus � −∞ ( ) → −(σ * − 2)/2 and x 2 ( ) → π when � ∞ → −(σ * − 2)/2 . Therefore, a critical undercooling of the dry GB may be defined as www.nature.com/scientificreports/ describing the fact that, even if σ * > 2 , the growth of the liquid layer along the dry GB does not occur when ∞ . Let us now consider the case corresponding to the other transition, i.e. the onset of melting for which x 2 = 0 and = 0 in Eqs. (10) and (11). Setting −∞ = 0, x 2 = 0 and x 1 = 2π in Eq. (15) and in the equilibrium profiles, we find a new critical value for ∞ : Fig. 2. For the typical values used previously and yielding L = 15 , it can be seen that � (2) ∞ > 0 . Thus, from our theoretical analysis of energy conservation during steady-state growth of the liquid layer, we find that premelting may take place even when the temperature to which the dry GB is brought lies above the melting temperature. We are not aware of any analysis in the literature pointing in that direction, probably because the dynamics in such a situation has rather scarcely been studied. In Refs. 35,37 , steady-state melting along a dry GB was investigated when σ * − 2 < 0 in an infinite system ( → ∞ ). Let us note that the steady-state velocity then diverges when σ * − 2 approaches 0.
The -dependence of the critical temperature � (2) ∞ implies that, for a given ∞ , the transition from premelting to melting may occur at This is what is shown in Figs. 7, 8, and 9, where we plot respectively � −∞ , x 2 /π and the degree of disorder of the liquid layer φ 3 (0) as a function of for a prescribed overheating ∞ = 0.023 , and for the same parameters as for Figs. 4, 5, and 6.
In Fig. 7, it can be seen that, in contrast to Fig. 4, −∞ does not converge to ∞ when → ∞ . Instead, −∞ reaches 0 at = c , and remains equal to 0 for larger . In Fig. 8, it can be seen that x 2 decreases from π to 0 when , smaller than the critical value, increases. Then, for > c , x 2 follows Here, in these plots, the fact that L is much larger than σ * − 2 explains why c is very similar for all three curves.
If now we fix , we may present −∞ as a function of ∞ . This is what is displayed in Fig. 10  To summarize the phenomenology associated with the two critical values of ∞ , we present in Fig. 11 a schematic kinetic phase diagram where, as a function of ∞ and , the equilibrium state after transformation is given. For � ∞ < � (1) ∞ , no transformation occurs. For � (1) ∞ < � ∞ < � (2) ∞ , the dry GB transforms into a premelted GB with x 2 > 0 and for � ∞ > � (2) ∞ melting occurs with x 2 < 0 . The interval of temperature for which the dry GB transforms into a pre-melted GB increases when decreases.
Let us note that, in the case where the disjoining potential's dependence on the layer's width is non-monotonous, for example yielding a thin-to-thick transition mentioned in the Introduction, corresponding complications of the scenarios presented above for −∞ , x 2 and φ 3 occur and turning points appear. In this case, the sign of the disjoining pressure, i.e. the derivative of the disjoining potential with respect to the layer's width, alternates, and −∞ is positive when the layer's width corresponds to an attractive part of the disjoining potential. We have nevertheless checked using a sharp-interface analytical approach using the disjoining potential proposed in Ref. 25 that � −∞ < � ∞ in any case. www.nature.com/scientificreports/ To conclude before presenting the simulations, the phase field model allows for a unified description including the microscopic physics at the atomic scale (represented by the scale of the interface width in the phase field model) responsible for premelting, and reproduces the macroscopic limit, for which the atomic distance formally vanishes, with the width of the liquid layer growing proportionally to when > c . Conversely, when is fixed, the width of the liquid layer decreases when ∞ , larger than � ∞ , decreases, and it reaches the atomic scale at � ∞ ∼ � (2) ∞ . As schematically represented in Fig. 11, the smaller , the larger ∞ should be in order for the liquid layer to achieve a macroscopic width. The difference between ∞ and −∞ (which exactly equals � (2) ∞ when � ∞ = � (2) ∞ since −∞ = 0 in this case) increases when decreases (with roughly a 1/ dependence). For the example displayed in Fig. 10 with = 4π , those differences are of order 0.1, corresponding to variations of temperature of order 20 K. These huge temperature differences are reached because is extremely small (approximately few tens of atoms). If we identify with the grain size of a polycrystalline structure, such a small is obviously unrealistic even for structures that are produced at extremely large cooling rates. For grains about few µ m in size, typically seen in laser-based additive manufacturing, we expect the temperature differences to be of the order of 0.1 K.

Phase field simulations
In this section, simulations of the growth dynamics of the liquid layer along a dry GB are presented. In order to obtain the correct conservation of energy in the simulations, we set the grid spacing to 1/75 of the interface width η and simulated domains that are very elongated in the growth direction z. From Fig. 10, it can be seen that the temperature of the equilibrated region far behind the tip obtained from the simulations (the dots) agrees well with the one calculated analytically (the solid lines). The diffusion coefficient is set constant throughout  Figure 11. Kinetic phase diagram where we give, as a function of ∞ and , the nature of the transformation: the GB remains dry when � ∞ < � (1) ∞ , the dry GB transforms into a pre-melted GB with an atomically thin liquid layer for � (1) ∞ < � ∞ < � (2) ∞ , and the dry GB transforms into a macroscopic liquid-solid equilibrium with a liquid layer having a macroscopic width when � www.nature.com/scientificreports/ the whole simulation domain, i.e. the diffusivity is equal in the liquid and solid phase. In Fig. 12a, where actually only a portion of the simulation is displayed and where the z axis is such that the steady-state velocity is negative, we present the map for the degree of disorder φ 3 (x, z) and for the dimensionless temperature field �(x, z) = Lη[T(x, z) − T m ]/(4σ sl T m ) (for a simulation at ∞ = −0.0046875 , σ * = 2.25 , = 4π and otherwise the same parameters as for Fig. 10, on which we denote −∞ obtained from the simulation using dots). In the region z < 0 , the GB is dry with φ 3 = 0 , and = ∞ . At small positive z, we distinguish a diffuse region where φ 3 emerges, and, as z increases the liquid layer achieves progressively its equilibrium structure. In Fig. 12b, the profiles of φ 3 along the x direction x, which is perpendicular to growth direction z, are displayed for z values exhibited in Fig. 12a. It can be seen that only at position z ≈ 300 , the liquid layer reaches its equilibrium structure, with the φ 3 profile changing very slightly between z = 276.5 and z = 360.2 . This is the distance to the tip of the liquid layer needed for the temperature field to equilibrate. It is obvious that simulating such a process using an atomistic method such as phase-field-crystal or molecular dynamics would be computationally intractable. In Fig. 13, the z-dependence of the phase field φ 3 and the dimensionless temperature averaged over the x-direction, i.e. �φ 3 �(z) = x (dx/ )φ 3 (x, z) and ���(z) = x (dx/ )�(x, z) , are presented. It can be seen that starts to deviate from ∞ and heat fluxes establish as soon as φ 3 deviates from 0. The latter deviation is linear, illustrating the fact that, as can be seen in Fig. 12b, the bump that characterizes φ 3 in the neighborhood of z = 0 has a width that varies weakly with z, while its amplitude varies much stronger. The saturation of φ 3 and at large z indicates that, far behind the tip, heat fluxes vanish and the system equilibrates.

Structure of the triple junction.
An important feature of the phase field results concerns the structure close to the tip of the liquid layer. Here, we focus only on the simulation corresponding to Fig. 12. In the macroscopic approach, an equilibrium in the form of Young's law at the triple junction between grain 1, grain 2 and the liquid layer is not possible for σ * − 2 > 0 , and the liquid should wet the GB with a vanishing macroscopic contact angle. Obviously in our case this picture is not appropriate since we study the propagation of the transition region that connects spatially the dry GB to the wet GB. We therefore focus on the behavior of the phase fields, www.nature.com/scientificreports/ shown for φ 3 in Fig. 12, on the scale of the interface width. We explore the structure close to the tip in a way similar to what is presented in Fig. 9 in Ref. 38 . We locate the two isolines φ 1 = φ 3 and φ 2 = φ 3 , and the intersection of these two isolines, at x = 0, z = z 0 , yields the position of the triple junction where φ 1 = φ 2 = φ 3 = 1/3 . The microscopic contact angle θ 0 is then defined by the angle between the isolines and the vertical axis at the triple junction. We illustrate this procedure in Fig. 14 and find θ 0 ≃ 2 • . In Ref. 35 , a generalization of Young's law for a triple junction at melting temperature is proposed in order to include effects of structural forces in the case σ * < 2 . The microscopic contact angle θ 0 relates the macroscopic contact angle θ ∞ and σ * through θ 2 ∞ − θ 2 0 = −(σ * − 2) . In the case σ * > 2 , their formula may actually be read by setting θ ∞ = 0 (according to classical Young's law), and one finds a finite microscopic contact angle θ 0 = √ σ * − 2 and a well-defined triple junction. We see that our simulation does not reproduce this result with a θ 0 much smaller than √ σ * − 2 . Several arguments may be invoked in order to explain this discrepancy. First, the analysis in Ref. 35 holds for a triple junction at melting temperature, while, in our simulation, the temperature field is highly inhomogeneous on the scale of the interface width and the temperature lies well below T m where the triple junction is defined. Second, their analysis is performed within a small slope approximation that requires σ * − 2 ≪ 1 , which is not the case in our simulation. Finally, as mentioned in Ref. 35 , kinetic effects at the triple junction may take place yielding a deviation from Young's law for the macroscopic contact angle. All these effects are interesting but their investigation is beyond the scope of this paper.  Figure 13. Dependence in the z-direction of the phase field φ 3 and the dimensionless temperature averaged over the x-direction (see text) for the simulation presented in Fig. 12.

Figure 14.
Isolines φ 1 = φ 3 and φ 2 = φ 3 in the neighborhood of the triple junction, located at ( x = 0, z = z 0 ) and corresponding to φ 1 = φ 2 = φ 3 = 1/3 . The contact angle θ 0 , given by the slope of the isolines at z 0 , is close to 2 • . www.nature.com/scientificreports/ We note also that the procedure outlined above and illustrated in Fig. 14 is enabled by the fact that φ 3 > 1/3 in the equilibrated region far behind the tip, as can be seen in Fig. 12. Thus, it is clear that when φ 3 remains smaller than 1/3, as for example for small enough , such a procedure may not be adopted. This case lies also beyond the scope of this paper.
Velocity dependence on surface energies. Let us now focus on the dependence of the tip growth velocity on the surface energies. For this purpose, we define the dimensionless velocity where V is the steady-state velocity and D is the heat diffusion coefficient. The normalization of the dimensionless velocity by 1/S comes from the following arguments. Our definition of in Eq. (8) differs from the usual definition of the dimensionless driving force for heat diffusion controlled processes in sharp interface models � sh = c p (T − T m )/L = S�/L . Since the growth velocity typically scales with a positive power of sh , it vanishes when c p /L → 0 , this fact being the reason why we divide our velocity by S. As we will see in the following, the definition of Ṽ given in Eq. (18) provides values of order unity. We set the channel width to = 4π , and we perform three simulations for each dependence, i.e. on σ * and on the dimensionless solid-liquid interface energy L −1 = σ sl /(Lη) . In Fig. 15, we plot the growth velocity as a function of σ * − 2 = 0.1, 0.2 and 0.25 for L = 15 and ∞ = 0 (we recall that, in a channel, ∞ = 0 does not represent a special point as illustrated by Fig. 10 for example). In Fig. 16, we plot velocity as a function of L −1 = 11/300, 1/15 and 2/15 for σ * = 2.25 and ∞ = −0.0046875 . We see that for both curves, the data is well fitted using a straight line passing through the origin. This suggests that V ∝ (σ * − 2)σ sl = σ gb − 2σ sl . Thus, it seems that in the premelting problem, the excess of interface energy attributable to the unfavorable dryness of the GB is the driving force. This kind of dependence of steady velocity on the interface energy is rather unusual for a growth process in a solid-liquid system. Indeed, in dendritic growth for example, the velocity is inversely proportional to the solid-liquid energy, www.nature.com/scientificreports/ the latter hindering the growth owing to the energetic penalization of interface curvature. Closer to our current interest, the steady-state velocity scales with the inverse of the solid-liquid interface energy also in the theories developed in Refs. 35,37 for melting along a GB. In the second part of the following discussion, we comment on these theories and on their link with our premelting problem. Let us comment on the values of the steady-state velocity that we obtain in phase field simulations. Using a typical value D = 10 6 m 2 /s for the heat diffusivity and the values for η and S used previously, i.e. η = 7.5 nm and S = 4c P T m /L = 40 , we find velocities of order 10 3 m s −1 . Such a value is questionable because it is in the order of the speed of sound. In our simulations, we have used a mobility M ij (see "Methods") providing equilibrium boundary conditions at the interface within the thin interface limit of the phase-field model 39 . Such a choice corresponds to an infinite interface mobility within the sharp interface limit for conditions close to equilibrium. Then, for sufficiently small velocities, the interface mobility does not play any role, neither in the pattern nor in the velocity selection mechanisms. Here, we have found that reducing the phase field mobility M ij does influence the growth. More precisely, reducing the phase field mobility M ij by a factor 30, approximately decreases the steady-state velocity by the same factor 30. This indicates that the kinetic does not follow the diffusion-controlled regime for which interface mobility does not play a role. On the other hand, for a purely kinetically-controlled regime, one should obtain a homogeneous temperature field, which here is not the case either (see Fig. 12). Thus, it is likely that the growth kinetics of pre-melted GBs falls into a mixed-mode regime where the interface mobility is a material related quantity that influences the transformation speed. This is known for other fast transformations, such as rapid crystallization in phase change materials 40 . The results may also be applied to premelting in an alloy. In this case, the diffusion coefficient in the liquid layer is probably in the order of a bulk liquid diffusion coefficient, i.e. typically three orders of magnitude smaller than the heat diffusivity provided above, and the diffusion coefficient in the solid grains is typically six orders of magnitude smaller than this heat diffusivity. If one assumes that the diffusion in the solid grains mainly sets the magnitude of the liquid layer growth velocity, the latter then becomes of the order of a mm/s, which is now well below the speed of sound. Let us note moreover that when one applies the classical dendritic growth theory to the melting of an alloy, for which the diffusion coefficient is much larger in the growing phase than in the disappearing phase, the velocity does not scale as D S but as D 2 S /D L where D S ( D L ) is the diffusion coefficient in the solid (liquid) 41 . This suggests that premelting in an alloy might even take place at velocities smaller than a mm/s.

Simulation of a polycrystalline evolution
Now, we briefly investigate the effect of premelting on the evolution of a polycrystalline structure relevant to experiments. We simulated the evolution of a microstructure that, initially, consisted of ten grains and a small liquid domain that enables the premelting liquid layer to propagate along the dry GBs. The GB energies are all equal with σ * = 2.25 , and as previously, the system is insulated. The evolution of the system as the liquid premelts the dry GBs and as the pre-melted GBs evolve with time can be seen in the "Supplementary Material". Here, we show only the results at three different points in time. In Fig. 17, an intermediate state for which the liquid layer has not yet propagated along every dry GB is displayed. In the left panel, the color codes φ 3 , in the middle panel a color is assigned to each grain, and in the right panel the color codes the dimensionless temperature . We see that the temperature has fallen below T m in the neighborhood of the liquid layer, while it is still at the initial temperature = 0.0046875 above T m in the bulk of the grains.
In Fig. 18, the liquid layer has propagated along all GBs; there has been an evolution of the size and shape of the grains; and we observe the heat diffusion field from the center of the grains to the pre-melted GBs. At the triple junctions, relatively large liquid pockets have emerged, where φ 3 = 1 since the distance between solid/ liquid interfaces is much larger than the interface width. Nevertheless, the curvature of the solid/liquid interfaces allows a temperature below T m in these regions. Interestingly, the protrusion that is highlighted with an arrow in the center panel is already observable at a much earlier stage of the evolution in Fig. 17. This is counterintuitive if, naively, one expects that the protrusion simply enhances the interface energy of the system. Instead, one www.nature.com/scientificreports/ should note that the curvature of the solid/liquid interface allows a minimization of the disjoining potential in the thin liquid layers adjacent to the triple junction, and is thus energetically favorable. This is an effect of the disjoining potential that should be of fundamental importance for the grain evolution in systems that are prone to premelting, and this point is dicussed further next. As a result of the promotion of premelting against melting due to confinement (here provided by small grain sizes) described in the previous sections, we see in Fig. 19 that the temperature everywhere in the domain is below T m although = 0.0046875 initially, in accordance with our analysis of steady-state growth of a liquid layer and the associated energy conservation (i.e., Eq. 15). At the stage of the evolution presented in Fig. 19, some grains have increased in size, some have decreased in size, and three grains have vanished. The liquid pockets at the triple junctions have further increased in size, but we still identify the effect of the disjoining potential minimization through unusual shapes of the grains and curvatures of the solid/liquid interfaces. Let us note that such premelted polycrystalline structure was nicely evidenced experimentally (see Figs. 3 and 4 in Ref. 42 ). These experiments were conducted in pure Al (that has, as mentioned earlier, properties close to the ones that we use here), and the undercooling ≃ −0.01 in our simulation, that corresponds to (T − T m )/T m = 0.01 × 4/L ≃ 0.0027 , is in qualitative agreement with the observation of the authors in 42 that "no signs of melting were detected for temperatures up to 0.999 T m .
In addition, we have plotted in Fig. 20 the φ 3 and profiles along the vertical line visible in the left and right panels of Fig. 19. We see that the temperature gradient is much larger across the liquid layer than in the solid. This gradient of temperature across the liquid layer, mainly resulting from the opposite sign of the Gibbs-Thomson effect at the two solid/liquid interfaces, is responsible for the motion of the liquid layer in its normal direction, with a velocity that is easily estimated using the Stefan condition at the interfaces.
Let us now present the physical picture that underpins the phenomenology of the polycrystalline evolution discussed above. First, note that a polycrystalline equilibrium state exists for � < 0 , i.e. at a temperature below T m , as shown in Fig. 21 (see the caption for more details). As already mentioned, when the size of the liquid pocket at the triple junctions is large enough, the liquid possesses the characteristics of the bulk liquid ( φ 3 = 1 ) at least at its center.  www.nature.com/scientificreports/ Of course, by symmetry, if the dry GB energies are equal, the angles between the pre-melted GBs will be the same (i.e., 2π/3 ); if those energies are not equal, the angles will still be the same if the radius of the curvature of the solid/liquid interfaces bounding the liquid pocket is much larger than the width of the layer, i.e., when | | ≪ 1 . Then, in an out-of-equilibrium scenario similar to the one discussed in Fig. 19, if grains are large enough, the curvature of the pre-melted layer will be small enough to make the temperature variations across the liquid layer much smaller than | | . In that case, the three lenght scales (the width of the pre-melted liquid layers η , their radius of curvature R of the order of grain size, and the radius of curvature of the liquid pockets η/|�| ) will be fully separated; the triple junction will be close to the equilibrium presented in Fig. 21; and the normal velocity of the pre-melted liquid layers will be proportional to their curvature, corresponding precisely to a curvature-driven phenomenon as for the commonly accepted dynamics of dry GBs. However, here, the mobility of pre-melted GBs, i.e. the proportionality coefficient between velocity and curvature, will be given by the diffusion coefficient in the liquid and will be therefore much higher than of dry GBs. Thus, in the regime η/R ≪ |�| ≪ 1 , the pre-melted GB dynamics corresponds to a classical curvature-driven GB dynamics with an enhanced mobility. Therefore, further investigations should focus on determining the structure of the triple junction, i.e. the structure of the liquid pocket, in these out-of-equilibrium conditions, and especially the deviation from the perfect 2π/3 angles present at equilibrium. Further investigations should also address the question of the existence of out-of-equilibrium steady-state situations.

Discussion
For our premelting problem, the phase field model is used in a way at odds with the usual philosophy of phase field theories, for example describing the pattern formation during solidification. Within the latter, the phase fields are functions that indicate which state the system presents locally. This state belongs to a discrete set, corresponding to the different stable or metastable phases under consideration. To each state corresponds a value of the phase fields, and the locus of their variations defines the interfaces. The interfaces' width η , much smaller than the length scale of the pattern and the diffusion length in the case of diffusion-controlled processes, is considered as a purely numerical parameter upon which no outcome of the simulation should depend. In this respect, the bulk of the different phases and the interfaces between them are well separated and identified. The interface width η should however be handled with care owing to its influence on the interface boundary conditions that the phase field model reproduces. Indeed, when the phase field model is recast into its corresponding sharp interface model using the so-called 'thin-interface limit' 39 , the interface boundary conditions depend on the equilibrium profiles of the phase fields and thus on η . The thin-interface limit is relevant at low enough driving forces, corresponding to linear kinetics obeying Onsager out-of-equilibrium thermodynamics 43 . Then, one is able to compare quantitatively the phase field results and the sharp-interface results, for example the growth velocity and the interfaces' shape.  www.nature.com/scientificreports/ In opposition, for our premelting problem, the characteristic length scale of the liquid layer is η , and therefore no separation of length scales holds. The bulk of the liquid phase and the solid-liquid interfaces cannot be identified separately, and the liquid layer corresponds to a multi-phase region where at least two phase fields do not vanish. Moreover a physical meaning is attributed to the profile of the liquid phase field φ 3 . For example, its maximum, that varies continuously with temperature according to Eq. (9), is denoted as the ' degree of disorder' .
The premelting dynamics that we have simulated is, to a large extent, governed by the coupling of the phase fields in the multi-phase region and a special care should thus be taken in order to discretize the phase fields' variations faithfully. As mentioned before, we had to use a grid step 75 times smaller than η in order to reach a good precision on the conservation of energy. This is a tremendous increase in the number of grid points used to discretize the interface compared to usual studies, for example in solidification. We are not aware of any example for which such a resolution was needed. This necessity restricts the size of the domains that we were able to simulate. As an illustration, since we had to use very elongated simulation boxes in the z-direction in order to quantitatively reproduce the analytical results for −∞ , we had to choose relatively small values of , i.e. 4π in the simulations presented in this article.
As we have seen previously, using a small promotes premelting against melting, which occurs only when a certain overheating of the dry GB is exceeded, i.e. for 0 < � (2) ∞ < � ∞ . In Refs. 35,37 , melting along a dry GB was studied using a sharp-interface approach. A crucial ingredient of the theory is the contact angle at the triple junction. In 37 , only the macroscopic contact angle θ ∞ provided by Young's law is considered. As mentioned before, it vanishes in our case σ * > 2 , and then, the steady-state growth velocity diverges. In Ref. 35 , a disjoining potential with an exponentially decreasing magnitude and a microscopic contact angle θ 0 were introduced. As mentioned in the previous section, θ 0 is related to the macroscopic angle and σ * through θ 2 0 = θ 2 ∞ + σ * − 2 , yielding a finite θ 2 0 = σ * − 2 when θ ∞ = 0 . This regularization allows the velocity to become finite because, in the scaling laws inherited from the theory in 37 , θ ∞ is replaced by θ 0 in 35 .
In both theories, the velocity is inversely proportional to σ sl for a given contact angle. This result opposes our simulations' suggestion that the velocity is proportional to σ sl for a fixed σ * . This opposition may be apprehended when noting that the driving force for premelting is the reduction of interface energy parametrized by σ gb − 2σ sl = (σ * − 2)σ sl , while the driving force for melting is the reduction of bulk free energy due to the stabilization of the liquid phase above T m . It would be interesting to investigate the transition from an interfaceenergy-dominated driving force to a bulk-energy-dominated one with phase field simulations. The theories of melting in 35,37 that yield a velocity inversely proportional to σ sl are derived for an infinite system in the direction perpendicular to growth (this is numerically possible because the Green's function method that is used allows to eliminate the diffusion field and reduces to a search for the interface shape). As a consequence of energy conservation, the interface shape then assumes Ivantsov parabolic asymptotics far behind the tip 44 . The distance between the solid-liquid interfaces thus diverges, and the disjoining potential vanishes. Within our analysis of the phase field model, we have seen that, in a channel also, the disjoining potential vanishes in the equilibrated region far behind the tip when > c . However, setting only slightly beyond c does seem appropriate in order to reach a regime where the driving force is dominated by the reduction of bulk free energy. It seems indeed natural to suppose that this regime holds when the width of the liquid film is much larger than the spatial range on which structural forces act, i.e. when the liquid film is much larger than η . Then conservation of energy in the channel tells us that the width of the liquid film in the equilibrated region far behind the tip is approximately given by ∞ (especially when L is large). Thus, for a study of the premelting to melting transition under the close-toequilibrium conditions assumed in Refs. 35,37 ( ∞ ≪ 1 ), we need a simulation box such that ≫ η/� ∞ ≫ η . In view of the extremely fine discretization required for an accurate realization of the complex coupling of the phase fields in the multi-phase region mentioned above, we understand that the investigation using phase field simulations of the cross-over between the regimes of bulk-energy-dominated and interface-energy-dominated driving force is rather challenging. The difficulty comes, of course, from the multi-scale nature of the problem. Of course, the multi-scale nature of the problem represents also a challenge numerically for the simulation of polycrystalline evolutions.

Conclusion
A multi-phase field model with obstacle potentials was used to study, for the first time, the dynamics of grain boundary premelting. In the model, the disjoining potential describes structural forces on the scale of the interface width (on the atomic scale in a sharp-interface approach) and decreases monotonically with the distance from the solid-liquid interfaces. The model was used to simulate the steady-state growth of a liquid layer along a dry GB in an insulated channel and the evolution of a pre-melted polycrystalline microstructure.
Our results show that a transition exists from a premelting transformation, which produces a pre-melted equilibrium state with a finite disjoining potential energy, to a melting transformation, which produces a macroscopic solid-liquid equilibrium with a vanishing disjoining potential energy. The results also reveal that, due to energy conservation, confinement, which is linked to grain size (or channel width), promotes premelting against melting, implying that a certain overheating of dry GBs above the bulk melting temperature should be exceeded for melting to occur; that overheating is found to be inversly proportional to grain size. Conversely, for a given overheating, melting occurs only at large enough grain sizes. Our computational results also show that premelting dynamics is governed by the reduction of the interface energy, with a velocity proportional to σ gb − 2σ sl .
We found that a polycrystalline equilibrium exists in which the triple junction takes the form of a liquid pocket with a macroscopic size. Realizing the presence of that equilibrium allowed us to gain novel insights into the evolution of pre-melted polycrystalline microstructures. That evolution consists of two stages: in the first stage, the liquid premelts the dry GBs, and then, in the second stage, the fully pre-melted GBs evolve. Due to the presence of the polycrystalline equilibrium, the disjoining potential plays a crucial role not only during the first Scientific Reports | (2020) 10:21074 | https://doi.org/10.1038/s41598-020-77863-9 www.nature.com/scientificreports/ stage, but also during the second one. For example, it allows the grains to have morphologies that, in the absence of the disjoining potential, would be energetically unfavorable. In addition, we find that, again due to the presence of the polycrystalline equilibrium, if grains within a premelted microstructure are large enough, in weakly out-of-equilibrium conditions, the dynamics may be recast into the well-known curvature-driven dynamic of dry GBs, with a mobility that is enhanced by a factor proportional to the ratio of liquid to solid diffusion coefficients. An interesting future study is investigating melting along a dry GB when σ gb − 2σ sl > 0 , and especially the transition between the regime driven by the reduction of interface energy, as in our simulations, and the one driven by the reduction of bulk free energy. Another interesting future study concerns the influence of the disjoining potential on the contact angles at the triple junction where the dry GB meets the two solid-liquid interfaces. Concerning the evolution of pre-melted polycrystalline microstructures, the perspective is to study the phenomenology arising from, as we have shown, the existence of three relevant length scales, i.e. the width of the liquid layer, the size of the liquid pocket at the pre-melted triple junction and the grain size.

Methods
The governing equations of the model can be outlined as follows. The total free energy reads where the local free energy density that depends on the spatial distribution of the three phase fields and the temperature reads where σ * = σ gb /σ sl tunes the tendency to premelting that occurs when σ * > 2.
In a dimensionless form and expressing lengths in units of η/(2π) , we have with The phase fields evolution equations are given, when the three phases are present (i.e. φ 1 φ 2 φ 3 = 0 ), by with (i, j, k) = (1, 2, 3), (2, 3, 1) or (3, 1, 2). When one of the phase field vanishes, for example φ k = 0 , then The functional derivatives are given by On the other hand, the temperature equation reads where D is the diffusivity, taken constant throughout the whole simulation domain, and with L = Lη/σ sl and S = 4c p T m /L.

Data availability
The data that support the findings of this study are available on request from the corresponding author [GB].