# Fileset

[Thermodynamically consistent multi-phase-field model incorporating passive external domains.pdf](https://mdr.nims.go.jp/filesets/8da09753-86a2-4c42-b333-3bdc42dca2e9/download)

## Creator

[Akimitsu Ishii](https://orcid.org/0000-0002-9261-4047), Yusuke Matsuoka, [Toshiyuki Koyama](https://orcid.org/0000-0001-7424-4858)

## Rights

[Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/)

## Other metadata

[Thermodynamically consistent multi-phase-field model incorporating passive external domains](https://mdr.nims.go.jp/datasets/e87270f7-0c15-4d3d-96de-3260e0a26d28)

## Fulltext

TF_Template_Word_Windows_2016Science and Technology of Advanced Materials: MethodsISSN: 2766-0400 (Online) Journal homepage: www.tandfonline.com/journals/tstm20Thermodynamically consistent multi-phase-fieldmodel incorporating passive external domainsAkimitsu Ishii, Yusuke Matsuoka & Toshiyuki KoyamaTo cite this article: Akimitsu Ishii, Yusuke Matsuoka & Toshiyuki Koyama (31 Mar2026): Thermodynamically consistent multi-phase-field model incorporating passiveexternal domains, Science and Technology of Advanced Materials: Methods, DOI:10.1080/27660400.2026.2653418To link to this article:  https://doi.org/10.1080/27660400.2026.2653418© 2026 The Author(s). Published by NationalInstitute for Materials Science in partnershipwith Taylor & Francis GroupView supplementary material Accepted author version posted online: 31Mar 2026.Submit your article to this journal View related articles View Crossmark dataFull Terms & Conditions of access and use can be found athttps://www.tandfonline.com/action/journalInformation?journalCode=tstm20https://www.tandfonline.com/journals/tstm20?src=pdfhttps://www.tandfonline.com/action/showCitFormats?doi=10.1080/27660400.2026.2653418https://doi.org/10.1080/27660400.2026.2653418https://www.tandfonline.com/doi/suppl/10.1080/27660400.2026.2653418https://www.tandfonline.com/doi/suppl/10.1080/27660400.2026.2653418https://www.tandfonline.com/action/authorSubmission?journalCode=tstm20&show=instructions&src=pdfhttps://www.tandfonline.com/action/authorSubmission?journalCode=tstm20&show=instructions&src=pdfhttps://www.tandfonline.com/doi/mlt/10.1080/27660400.2026.2653418?src=pdfhttps://www.tandfonline.com/doi/mlt/10.1080/27660400.2026.2653418?src=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1080/27660400.2026.2653418&domain=pdf&date_stamp=31%20Mar%202026http://crossmark.crossref.org/dialog/?doi=10.1080/27660400.2026.2653418&domain=pdf&date_stamp=31%20Mar%202026https://www.tandfonline.com/action/journalInformation?journalCode=tstm201 Publisher: Taylor & Francis & The Author(s). Published by National Institute for Materials Science in partnership with Taylor & Francis Group Journal: Science and Technology of Advanced Materials: Methods DOI: 10.1080/27660400.2026.2653418 Thermodynamically consistent multi-phase-field model incorporating passive external domains Akimitsu Ishiia*, Yusuke Matsuokaa and Toshiyuki Koyamaa aResearch Center for Structural Materials, National Institute for Materials Science *Corresponding author: ISHII.Akimitsu@nims.go.jphttps://crossmark.crossref.org/dialog/?doi=10.1080/27660400.2026.2653418&domain=pdf 2  Thermodynamically consistent multi-phase-field model incorporating passive external domains Phase-field (PF) modeling of multiphase and multicomponent material processes where material phases and their external domains coexist, such as sintering, lacks a thermodynamically consistent formulation. To address this issue, this paper introduces a novel multi-phase-field (MPF) model with passive external domain (PED), which is defined as an external region that remains inert and does not contribute to microstructural evolution in the material system. An additional variable associated with concentration fields is incorporated to distinguish between the material and PED without thermodynamic inconsistencies. The model is validated through 2D single-phase simulations, which highlight its ability to capture surface triple junction behavior, consistent with Young’s equation. Its capability is further demonstrated through 3D simulations of a three-phase, ternary system coexisting with a PED. These simulations correctly predict the evolution of the system toward thermodynamic equilibrium, including phase transformations and solute partitioning. Overall, the MPF-PED model provides a robust and physically consistent framework, overcoming the limitations of conventional models in handling external domain to facilitate the realistic simulation of complex material processes. Keywords: Phase-field modeling, Microstructure, Multicomponent, Diffusion, CALPHAD 1: Introduction The phase-field (PF) method is a powerful numerical technique for simulating the evolution of complex microstructures, particularly those involving moving interfaces. As this method is fundamentally based on minimizing the total free energy of a system, it provides a unified framework for analyzing various physical phenomena. Consequently, the PF method has been widely applied to predict microstructural evolution in diverse materials, leading to the development of numerous PF models and applications. A significant advantage of the PF method is its ability to couple  3  thermodynamically consistently with concentration fields [1], which is crucial for materials design. Integration with the Calculation of Phase Diagrams (CALPHAD) method [2] further enables the prediction of microstructure evolution in realistic multiphase and multicomponent systems. However, applying PF simulations to practical materials design often requires accounting for an external domain, such as a vacuum or gas phase outside the material itself. This consideration is essential in processes such as sintering [3, 4], additive manufacturing [5–7], physical vapor deposition [8, 9], and electromigration [10–12]. In simple systems—particularly those involving a single-phase, for which analyzing concentration fields is intrinsically superfluous—a density field [3–5] or concentration field (virtual, in some cases) [13–16] is typically used to distinguish the material domain from the exterior. The variable representing these fields is usually treated as a conserved quantity and solved using the Cahn–Hilliard equation [17], ensuring mass conservation of the material. However, for multiphase–multicomponent systems, this problem remains unresolved owing to critical theoretical and computational limitations. This study focuses on external domains that remain inert and do not contribute to microstructural evolution in the material system, referring to such domains as passive external domains (PEDs). Conventional methodologies for incorporating PED in multiphase–multicomponent systems can be broadly categorized into two approaches: (i) Treating the external domain as a gas phase [18–21]: This approach is conceptually sound but practically flawed. It often relies on virtual Gibbs energies owing to the lack of CALPHAD databases, necessitating physically arbitrary adjustments to match established phase diagrams under varying process conditions. Furthermore, the significant disparity in properties (e.g., molar volume and diffusivity) between gas and condensed phases frequently leads to computational instability. (ii) Treating the external  4  domain as a vacuum [22]: While this approach avoids gas-phase complexities, it fails to provide a thermodynamically consistent framework for concentration fields because a vacuum is inherently immiscible with material phases and concentration undefined within it. Moreover, in conventional PF frameworks, the material–vacuum interface would need to be treated as if it were in equilibrium, a notion that is physically meaningless. Therefore, if the external domain of the material plays no significant role in the overall state change of the system and can be treated as a PED, its exclusion from the calculation is a reasonable and numerically robust approach. To overcome these challenges, this paper introduces a novel PF model capable of thermodynamically consistent treatment of PED. The core strategy involves introducing a new variable field that defines the domain where concentration fields are valid as independent variables. A pivotal advantage of this formulation is that the Gibbs energy of the PED is fundamentally excluded from the theoretical framework, thereby eliminating the need for nonphysical parameters while maintaining thermodynamic consistency at the interface. As the proposed model builds on the multi-phase-field (MPF) method [23, 24], it is termed the multi-phase-field with PED (MPF-PED) model. The MPF model is a promising approach for simulating polycrystalline phenomena such as grain growth [25–27], and various related methodologies have been established for analyzing the evolution of microstructure and concentration in multiphase–multicomponent systems [28–31]. Recent advancements have extended this capability to systems containing more than 20 components [32, 33]. Another significant benefit is the availability of well-developed algorithms for computational efficiency [34], enabling large-scale simulations within acceptable timescales [35, 36]. Building upon this framework, the MPF-PED model inherits these advantages while addressing the critical limitation associated with external domains.  5  This paper first presents the formulation of the MPF-PED model. The model is validated by examining the morphological and concentration evolution at the interface between a material phase and a PED in two dimensional (2D) single-phase simulations. Subsequently, the capability of the MPF-PED model is demonstrated through a three-dimensional (3D) simulation of polycrystalline grain growth in a three-phase, ternary system coexisting with PEDs. The MPF-PED model represents the first approach to overcome the limitations imposed by PEDs in PF simulations in a physical consistent manner.  2: Model description 2.1. MPF-PED model A microstructure including a PED is represented by two types of PF variables [37, 38]:  and  ( = 1, 2, … N). Figure 1 illustrates the profiles of these variables for a system containing two material phases and PEDs. The variable  distinguishes between the material domain ( = 1) and PED ( = 0), whereas  represents the probability of existence of solid or liquid phases. Although  is intrinsically undefined within the PED, the present model assumes that the PED retains  the value corresponding to the nearest material phase. Therefore, at the surface region,  changes smoothly from 0 to 1, whereas  remains constant. Conversely, in the internal interface regions (i.e., grain and phase boundaries) of the material,  = 1 and  transition smoothly. Hence, the domain of the ath phase is represented by the product . Additionally, ci (i = 1, 2, … Nc) denotes the molar fraction field variable of substitutional components. Both  and ci  6  are conserved variables, whereas  is a non-conserved variable.  and ci satisfy the following constraints: ∑           (1) and ∑            (2) The molar fraction ci is defined as    ∑            (3) where    represents the local molar fraction of the ith component in the th phase. To describe the concentration within the material domain, an effective concentration variable,  ̃  is defined as  ̃       (4) From Eqs. (2) and (4), the following relationship holds for  ̃  and :   ∑  ̃        (5) The physical interpretation of  ̃  can be clarified using the definitions of  and ci. Let n0 represent the total number of moles within a reference volume when it is fully occupied by a material, and n represent the actual total number of moles locally present within the same volume. Using n0 and n, the variable  can be expressed as  = n / n0, indicating the probability of existence of the material. The local molar fraction ci is defined as ci = ni / n, where ni is the number of moles of the ith component within the reference volume. Consequently,  ̃  = ci represents ni normalized by the reference molar density (n0). This formulation ensures that  ̃  correctly becomes zero in the PED (where n = 0), circumventing the physical ambiguity of defining ci where no material exists.  7  The total free energy of the system, G, is defined as   ∫ [                ]     (6) where f and f are the energy densities related to the distributions of  and , respectively; and fchem is the chemical free energy density. The interpolation function P() is defined as P() = 2(3 − 2). Detailed justification for this selection is provided in Supplementary Material S1. The energy density f is defined as the sum of the double-well potential and gradient energy density for :                         (7) where W and  denote the height of the energy barrier and gradient coefficient, respectively. Although the surface energy s can be expressed as a function of  to account for phase-dependent energies [39], it is assumed constant in this study for simplicity, independent of the material phase.  Accordingly, W = 12s /  and  = 3s / 2, where  is the thickness of the surface region [37]. The energy densities f and fchem adopt the forms commonly used in the MPF model:    ∑ ∑ (                   )            (8) and       ∑          (9) where W and  denote the height of the double-obstacle potential and gradient coefficient at the interface between the ath and bth phases, respectively. These are defined as W = 4 /  and  = 8 / 2, where  is the width of the internal interface region and  is the interfacial energy between the th and th phases. f denotes the chemical free energy density of the th phase, typically obtained from a CALPHAD database.  8  The evolution of  is directly coupled to that of  ̃ . Its time evolution equation is obtained via the time derivative of Eq. (5):         (∑  ̃      )  ∑  ̃          (10) Therefore, in the proposed MPF-PED model, the evolution of  is solved implicitly: It is determined from the time evolution of  ̃  for i = 1 to Nc. The time evolution equation of  ̃  is expressed using the following diffusion equation:   ̃      (       ̃ )   (11) where Mi is the diffusion mobility of the ith component. To account for different diffusion pathways, Mi is defined as                   (12) where    ,    , and      . are the volume diffusion, surface diffusion, and internal interfacial diffusion mobilities, respectively. Expressing each mobility as a function of the variables  and , bulk, surface, and interface locations and the phase types with their respective mobilities [40–43].         ∑                       ∑                ∑ ∑                 (13) where      is the volume diffusion mobility in the th phase, which can be determined using the atomic mobility data in the CALPHAD database [44];      denotes the surface diffusion mobility in the th phase; and     is the diffusion mobility at the internal interface between the th and  th phases. Q() is the interpolation function introduced to ensure numerical stability. In this study, Q() is defined as Q() = 3(10 − 15 + 62) as in the previous studies [3,5,40,41]. The weighting coefficients for surface and internal interfacial diffusion mobilities are normalized to achieve unity at the center of the surface region ( = 0.5) and at the centers of the internal interface region  9  ( =   = 0.5), respectively. Note that    = 0 when  = 0, which ensures that the mass flux is zero in the PED. This definition of mobility allows for a slight residual  in the wake of the interface migration. The preliminary study confirms that this residue is less than 1% of the value within the material phase. Since diffusion is suppressed in this low- region, the impact of this artifact on mass transport and mass conservation is negligible in most cases. However, we cannot rule out the possibility that this artifact may affect simulation results, particularly in scenarios involving extensive interface migration or when the interface width is set larger than conventional choices. Under such conditions, the incorporation of an additional term, such as an anti-trapping current [45], would be required. Considering that the third term on the right-hand side of Eq. (13) and that the relation ∑ = 1 holds, Eq. (11) can be rewritten as follows:   ̃      ∑ (          ̃ )      (14) with                                  ∑ ∑                 (15) In the MPF-PED model, the local equilibrium condition (i.e., the parallel tangent law) is assumed at the interface. If     can be treated as a dependent variable, the following condition of equal diffusion potentials holds [1]:  ̃                                      (16) However, this condition is not applicable to Eq. (11), as it independently solves for Nc concentrations. The functional derivative of the total free energy with respect to  ̃  is given by (see Supplementary Material S2 for derivation)  10      ̃  (                )            ∑     (   ∑                  ∑ ∑                       )  (17) The third term on the right-hand side of Eq. (17) reduces to the chemical potential of the ith component. The partial derivative ∂   / ∂ci can be derived analytically, as detailed in the Supplementary Material S3. In the MPF-PED model, the time evolution is solved for all Nc effective concentration variables scaled by . Simultaneously, their sum— itself—is driven toward 0 or 1 by a double-well potential acting as a soft constraint. Consequently, the correlations among all components are inherently incorporated through Eq. (11). The time evolution of  is described by [23, 24]               ∑     [          ∑ {(       )     (         )    }             √    {      ∑   (       )       }]                 (18) where Nlocal is the local number of phases with  > 0.   is the PF mobility, defined as   = (2/8) M, where M is the internal interfacial mobility between the th and th phases. For solving this equation,     can be treated as a dependent variable. The change in    over one time step,    , is calculated using the recently proposed direct CALPHAD coupling approach [32, 33]. By defining ’ and ’ as the values after one time step, the following relationship holds based on Eqs. (3) and (4):  ∑           ̃    ∑                     (19) where k is the ratio between     and    , given by [32, 33]  11      (               )                                             (20) Substituting Eq. (20) into Eq. (19) yields the following equation for the material domain ( > 0):         ∑      (               )                      ̃  ∑                   ∑      (                       ) (21) Then, Eq. (21) is solved for     (i ≤ Nc − 1). The dependent concentration change,     , is determined from the constraint in Eq. (2). By adopting this approach, the concentration change satisfying the equal diffusion potential condition is obtained without iterative calculations, allowing solute partitioning occurs at interface between material phases. It should be noted that, at the surface described in the MPF-PED model, only  changes while  remains constant. Consequently,     becomes zero; therefore, solute partitioning at the surface does not occur in principle.  2.2. Numerical Implementation Two numerical techniques were used to ensure stable simulations. The first method was aimed at preventing numerical overshooting/undershooting of . Although  is defined in the range of 0 to 1,  can overshoot or undershoot the defined range due to various numerical factors. A representative example arises from the relationship between  and  ̃ ; strong Gibbs–Thomson effects may cause numerical artifacts resulting in  > 1. While logarithmic or double-obstacle potentials could be employed to impose stricter constraints and prevent such bound violations, they are prone to causing numerical  12  instability when solving the Cahn–Hilliard equation. Therefore, to mitigate this issue, a conditional outer penalty was applied via the double-well potential coefficient in Eq. (7). While the coefficient retains the value W when 0 ≤  ≤ 1, it is typically set to a sufficiently large value in the numerical artifact region where  < 0 and  > 1. As this outer penalty term does not correspond to a physical material property, its magnitude can be chosen relatively freely as long as it fulfills the necessary requirements—namely, effectively suppressing bound violations without inducing numerical instability. In this study, the outer penalty coefficient     was set as     = 100W through a preliminary examination (see Supplementary Material S4). The second method was designed to suitably handle the PF variables, , within the surface region. In the MPF-PED model,  lacks physical meaning in the PED. Therefore, values of  from the PED must be excluded when calculating their spatial derivatives. As the equilibrium profile of  at the surface follows a hyperbolic tangent function in practice,  values were considered invalid in regions where  falls below a small threshold . Therefore, the zero Neumann boundary condition was applied to each  at the boundary defined by  = . Based on preliminary investigations,  was set as 0.0001. Notably, this boundary condition is consistent with the assumption that the PED inherits the phase state of the adjacent material, which implies a vanishing spatial derivative of  across the threshold . For spatial points where  increases and exceeds the threshold  due to interface migration, the values of  are initialized. Specifically,  at such points is determined by averaging the values from the immediate neighborhood where  > . A normalization step is then applied to these points to ensure ∑ = 1. Consequently, the summation constraint of Eq. (1) is strictly satisfied throughout the effective calculation domain.  13  The time evolution equations for a and  ̃  were discretized using the finite difference method (FDM). The first-order Euler FDM was used for time integration, and the second-order central FDM was used for spatial discretization.  3: Simulation conditions To validate the MPF-PED model, simulations were performed on a hypothetical ternary material system composed of elements A, B, and C. The chemical free energy density f for each phase was defined using the regular solution approximation:       (∑[                    ]      ∑ ∑                    ) (22) where Vm is the molar volume, R is the gas constant, T is the absolute temperature, 0   is the molar Gibbs energy of the ith solute component in the ath phase, and ij denotes the interaction parameters between the ith and jth components. The hypothetical system was designed to have three stable phases (X, Y, and Z) arranged symmetrically in the phase diagram at 1273 K (Figure 2). The equilibrium compositions of phases X, Y, and Z corresponded to molar fractions [A, B, C] of [0.8, 0.1, 0.1], [0.1, 0.8, 0.1], and [0.1, 0.1, 0.8], respectively. These compositions were obtained by setting the molar Gibbs energies to 0   = 3595.28 × (-1)   J/mol, where i is the Kronecker delta, and the interaction parameters to ij = 21170 J/mol (for i ≠ j). The simulation parameters are summarized in Table 1, expressed in non-dimensional forms normalized by the finite-difference grid size  l, internal interfacial mobility M, and energy density, RT/Vm. The molar volume was assumed to be uniform for all phases, with a value of Vm = 5 × 10−6 m3/mol. The simulations were conducted under isothermal conditions at 1273 K.  14   4: Results and Discussion 4.1. Model validation In the MPF-PED model, the material surface is represented by the variable , a feature absent in conventional MPF models. Therefore, validation of the morphological and concentration evolution at surfaces is required. In particular, the dihedral angle at the triple junctions formed by the surface and an internal interface serves as a benchmark [46]; agreement with Young’s equation demonstrates the model’s consistency with analytical solutions for interfacial energy balance. Notably, if the entire computational domain remains constant at  = 1, the MPF-PED model becomes equivalent to the conventional MPF model. Therefore, validation of the microstructural evolution within the bulk material was omitted, as the model inherently reproduces the conventional MPF behavior. Instead, 2D simulations of a single-phase system were conducted for quantitative verification.  Figure 3(a) shows the initial conditions for the simulations. Within the computational domain of 256 × 128 computational grids, two circular grains of phase X were placed adjacent to each other. Here, to describe the distributions of phases and interfaces, a field parameter  was defined for visualization:  Φ  ∑            (23) Regions with  = 0 denote PEDs; those with  = 1 denote the grain interiors; and those with 0 <  < 1 denote the interfaces. The interior of each grain initially possessed the equilibrium concentration. Using these initial conditions, simulations were performed for 2 × 106 computational steps while varying the ratio of surface energy to grain  15  boundary energy, s / XX, from 0.2 to 4.0. According to Young’s relation, if s / XX < 0.5, the grain boundary should be completely wetted by the PED. Therefore, the condition s / XX = 0.2 was simulated to confirm that the initially contacted grains are correctly separated due to the energetic drive for wetting, consistent with the theoretical prediction. Figure 3(b) presents the results obtained for s / XX = 1.0. The contact area (neck) between the grains grew over time, driven by the balance of interfacial energies at the triple junction. Figure 3(c) shows the profiles of  and  ̃  along the white dashed line indicated in Figure 3(b). As expected for a system evolving from an equilibrium state, the material within the grown neck maintained the equilibrium concentration of phase X. Figure 3(d) illustrates the resulting triple junction morphologies and dihedral angles obtained from simulations with seven different s / XX ratios: 0.2, 0.75, 1.0, 1.5, 2.0, 3.0, and 4.0. The dihedral angle , defined as the angle between the two surfaces, was evaluated using the surface normal vectors near the triple junction:            (                            |         |) (24) where the ()left and ()right represent the spatial gradients of  at the left and right vicinities of the triple junction, respectively. In this study, the vectors ()left and ()right were programmatically defined and calculated as follows. Focusing on the upper triple junction in the two-particle model, the calculation for the right side was performed first. In the right half of the computational domain (x ≥ 129 l), a reference grid point was identified by finding the location that satisfies ∑( < 0.6 and  > 0.5, specifically selecting the point where  is closest to 0.5. Based on this reference point, the spatial gradient  was calculated at each grid point within a 3 × 3 region extending in the upper-right direction, and their average was taken. To obtain the gradient for the left side, the same operation was applied symmetrically in a mirrored manner to the left  16  half of the domain. Finally, these averaged spatial gradients were assigned to ()left and ()right in Eq. (24). These numerical solutions and the analytical solutions calculated using Young’s equation exhibited close agreement for 0.75 ≤s / XX ≤ 4.0. In the case of s / XX = 0.2, the two phase-X grains that were initially in contact separated over time. These results demonstrate that the MPF-PED model can accurately analyze the evolution of interfacial triple junctions satisfying the interfacial energy balance.  4.2. Multiphase and multicomponent grain growth with PEDs A key feature of the MPF-PED model is its ability to simulate microstructure evolution in systems where multiple phases coexist with PEDs, without defining the Gibbs free energy for the PED. To validate that the predicted phase transformations and concentration evolution are consistent with the phase diagram under these conditions, a 3D simulation of grain growth was performed, involving the coexistence of three material phases (X, Y, and Z) and a PED. Spherical grains of these phases were randomly distributed within a computational domain of (128 l)3 grids under periodic boundary conditions. The grain radii for phases X, Y, and Z were set as 20 l, 15 l, and 10 l, respectively, and the corresponding numbers of grains were 12, 28, and 96. These values were chosen to ensure that the initial volumes of all phases were nearly identical. The total volume of the initial grains occupied 60% of the computational domain. The interfacial energy ratio was set as s /  = 2.0. Two simulation cases were designed, with different initial phase compositions. In Case 1, all phases initially possessed their equilibrium compositions (Figure 2). In Case 2, phase X was set to a non-equilibrium  17  state with (  ̃A ,  ̃ ,  ̃ ) = (0.6, 0.2, 0.2), whereas phases Y and Z remained in equilibrium. Figure 4(a) show snapshots of the morphological evolution for Case 1. As discussed in the previous section, the grain shapes evolved, and the initial necks between grains grew, driven by interfacial energy. Additionally, grain boundary migration driven by the interface curvature was observed between grains of the same phase, resulting in changes in individual grain sizes, although this effect was not prominent in this simulation. As all phases were initially in equilibrium, no driving force for phase transformation acted at the phase interfaces. Therefore, the volume of each phase was expected to remain constant. To quantitatively verify this, the evolution of phase volume fractions was evaluated. Figure 5(a) shows the temporal variation of these volume fractions. Dashed lines in the figure indicate the equilibrium phase volume fractions calculated using the CALPHAD method based on the overall composition. Although slight deviations from the theoretical values were observed due to the initial adjustment driven by curvature effects and minor numerical errors in the initial configuration, the simulated volume fractions remained constant and closely matched the calculated equilibrium fractions. At 5 × 105 computational steps, the relative errors between the simulated volume fractions and the equilibrium values were 0.37%, −0.03%, −0.04% for phases X, Y, and Z, respectively. Figure 4(b) illustrates the morphological evolution for Case 2. The X phase, initially in an unstable non-equilibrium state, decreased in volume over time. Concurrently, the surrounding Y and Z phases grew. The evolution of the concentration field  ̃A within the material domain is depicted in Figure 4(c). The concentration  ̃A inside the X phase increased, approaching its equilibrium value. This concentration change was attributable to solute partitioning at the interfaces between the X phase and  18  other phases (Y and Z). By definition of the MPF-PED model, solute partitioning does not occur at the interface with PED in principle. However, elevated concentrations were observed at the surface of the X phase at 5 × 105 steps, attributable to rapid surface diffusion relative to bulk diffusion. By 5 × 105 steps, the interior of the X phase had approached its equilibrium concentration of  ̃A = 0.8. It should be noted that the current model assumes a constant surface energy s and does not account for the Gibbs adsorption effect, where surface energy varies with local concentration. While the observed concentration accumulation at the surface is a kinetic result of rapid surface diffusion rather than thermodynamic segregation, future extensions of the model would be required to study systems where the composition-dependence of surface energy plays a critical role. Figure 5(b) shows the evolution of the phase volume fractions for Case 2. Consistent with the morphological changes, the volume of the X phase decreased, whereas those of the Y and Z phases increased. The convergence of the volume fractions toward their equilibrium values by 5 × 105 steps is consistent with the observation that the concentration within the X phase reached equilibrium. Similarly, at 5 × 105 computational steps, the relative errors were 1.25%, −0.27%, and −0.50%, for phases X, Y, and Z, respectively. These results demonstrate that the proposed MPF-PED model can accurately analyze microstructural evolution in multiphase, multicomponent systems containing PEDs, while maintaining thermodynamic consistency.   19  5: Conclusions This study proposed and validated a novel MPF-PED model designed to appropriately treat PEDs in simulations of multiphase–multicomponent materials. The proposed approach enables thermodynamically consistent analyses by introducing an additional variable associated with the concentration fields into the conventional MPF framework, thereby effectively distinguishing the material domain from the PED. The model was first validated through 2D single-phase simulations, which demonstrated that the MPF-PED model accurately captures the interfacial energy balance at triple junctions involving the material surface, as well as the dynamics of surface-involved grain growth. Furthermore, the model’s capability was verified through 3D simulations of a complex three-phase, ternary system coexisting with a PED. The results showed that the model successfully captures the evolution toward thermodynamic equilibrium, accurately predicting both microstructural changes and solute partitioning. Despite its demonstrated capabilities, it should be noted that the current formulation is presently restricted to handling substitutional elements. Future development will be required to extend the framework to incorporate interstitial elements. Overall, the proposed MPF-PED model provides a physically consistent and robust framework for analyzing systems with external domains, thereby overcoming a critical limitation of conventional PF models. The MPF-PED model is therefore expected to be applicable to microstructure analyses in various material processes, such as sintering, additive manufacturing, and physical vapor deposition. As a specific direction for future work, we plan to apply this model to the detailed analysis of sintering phenomena, where the thermodynamically consistent treatment of evolving  20  PEDs is critical. This represents a significant advancement toward the quantitative and thermodynamically consistent simulation of these complex material processes.  Acknowledgements This study was supported by the Grant-in-Aid for Early-Career Scientists (JSPS KAKENHI, Japan) [JP24K17179] and the MEXT Program: Data Creation and Utilization-Type Material Research and Development Project (Digital Transformation Initiative Center for Magnetic Materials) [JPMXP1122715503].  References [1]  Kim SG, Kim WT, Suzuki T. Phase-field model for binary alloys. Phys Rev E. 1999;60(6):7186–7197. doi:10.1103/PhysRevE.60.7186. [2]  Ågren J. Calculation of phase diagrams: Calphad. Curr Opin Solid State Mater Sci. 1996;1(3):355–360. doi:10.1016/S1359-0286(96)80025-8. [3]  Wang YU. Computer modeling and simulation of solid-state sintering: A phase field approach. Acta Mater. 2006;54(4):953–961. doi:10.1016/j.actamat.2005.10.032. [4]  Xue M, Yi M. Phase-Field Simulation of Sintering Process: A Review. Comput Model Eng Sci. 2024;140(2):1165–1204. doi:10.32604/cmes.2024.049367. [5]  Yang Y, Ragnvaldsen O, Bai Y, et al. 3D non-isothermal phase-field simulation of microstructure evolution during selective laser sintering. npj Comput Mater. 2019;5(1):81. doi:10.1038/s41524-019-0219-7. [6]  Lu L-X, Sridhar N, Zhang Y-W. Phase field simulation of powder bed-based additive manufacturing. Acta Mater. 2018;144:801–809. doi:10.1016/j.actamat.2017.11.033. [7]  Keller T, Lindwall G, Ghosh S, et al. Application of finite element, phase-field, and CALPHAD-based methods to additive manufacturing of Ni-based superalloys. Acta Mater. 2017;139:244–253. doi:10.1016/j.actamat.2017.05.003. [8]  Dai R, Yang S, Zhang T, et al. High-Throughput Screening of Optimal Process Parameters for PVD TiN Coatings With Best Properties Through a Combination of 3-D Quantitative Phase-Field Simulation and Hierarchical Multi-Objective Optimization Strategy. Front Mater. 2022;9:924294. doi:10.3389/fmats.2022.924294.  21  [9]  Stewart JA, Spearot DE. Phase-field simulations of microstructure evolution during physical vapor deposition of single-phase thin films. Comput Mater Sci. 2017;131:170–177. doi:10.1016/j.commatsci.2017.01.034. [10]  Guo Y, Huang P. A multi-phase-field model of void crossing grain boundary under electromigration-induced anisotropic surface diffusion in interconnects. Eur J Mech A Solids. 2024;106:105305. doi:10.1016/j.euromechsol.2024.105305. [11]  Santoki J, Mukherjee A, Schneider D, et al. Effect of conductivity on the electromigration-induced morphological evolution of islands with high symmetries of surface diffusional anisotropy. J Appl Phys. 2021;129(2):025110. doi:10.1063/5.0033228. [12]  Ishii A, Yamanaka A. Multi-phase-field modelling of electromigration-induced void migration in interconnect lines having bamboo structures. Comput Mater Sci. 2020;184:109848. doi:10.1016/j.commatsci.2020.109848. [13]  Nakazawa A, Sakane S, Takaki T. Multi-phase-field modeling and high-performance computation for predicting material microstructure evolution during sintering. J Mater Res Technol. 2025;34:1803–1816. doi:10.1016/j.jmrt.2024.12.171. [14]  Greenquist I, Tonks MR, Aagesen LK, et al. Development of a microstructural grand potential-based sintering model. Comput Mater Sci. 2020;172:109288. doi:10.1016/j.commatsci.2019.109288. [15]  Hötzer J, Seiz M, Kellner M, et al. Phase-field simulation of solid state sintering. Acta Mater. 2019;164:184–195. doi:10.1016/j.actamat.2018.10.021. [16]  Chakraborty S, Kumar P, Choudhury A. Phase-field modeling of grain-boundary grooving and migration under electric current and thermal gradient. Acta Mater. 2018;153:377–390. doi:10.1016/j.actamat.2018.04.019. [17]  Cahn JW, Hilliard JE. Free Energy of a Nonuniform System. I. Interfacial Free Energy. J Chem Phys. 1958;28(2):258–267. doi:10.1063/1.1744102. [18]  Ishii A, Koyama T, Abe T, et al. Multi-phase-field modeling of sintering applicable to solid-state and liquid-phase sintering in multiphase and multicomponent systems. Mater Today Commun. 2024;40:110116. doi:10.1016/j.mtcomm.2024.110116. [19]  Villanueva W, Grönhagen K, Amberg G, et al. Multicomponent and multiphase simulation of liquid-phase sintering. Comput Mater Sci. 2009;47(2):512–520. doi:10.1016/j.commatsci.2009.09.018. [20]  Wang F, Nestler B. A phase-field study on the formation of the intermetallic Al2Au phase in the Al–Au system. Acta Mater. 2015;95:65–73. doi:10.1016/j.actamat.2015.05.002. [21]  Wang F, Reiter A, Kellner M, et al. Phase-field modeling of reactive wetting and growth of the intermetallic Al2Au phase in the Al-Au system. Acta Mater. 2018;146:106–118. doi:10.1016/j.actamat.2017.12.015.  22  [22]  Ishii A, Koyama T, Abe T, et al. Phase-field simulation coupled with a CALPHAD database for analyzing microstructural evolution during liquid-phase sintering of Nd–Fe–B magnets. J Alloy Compd. 2025;1025:180266. doi:10.1016/j.jallcom.2025.180266. [23]  Steinbach I, Pezzolla F, Nestler B, et al. A phase field concept for multiphase systems. Physica D. 1996;94(3):135–147. doi:10.1016/0167-2789(95)00298-7. [24]  Steinbach I, Pezzolla F. A generalized field method for multiphase transformations using interface fields. Physica D. 1999;134(4):385–393. doi:10.1016/S0167-2789(99)00129-3. [25]  Miyoshi E, Takaki T, Ohno M, et al. Accuracy Evaluation of Phase-field Models for Grain Growth Simulation with Anisotropic Grain Boundary Properties. ISIJ Int. 2020;60(1):160–167. doi:10.2355/isijinternational.ISIJINT-2019-305. [26]  Miyoshi E, Takaki T. Multi-phase-field study of the effects of anisotropic grain-boundary properties on polycrystalline grain growth. J Cryst Growth. 2017;474:160–165. doi:10.1016/j.jcrysgro.2016.11.097. [27]  Shahnooshi E, Jamshidian M, Jafari M, et al. Phase field modeling of stressed grain growth: Effect of inclination and misorientation dependence of grain boundary energy. J Cryst Growth. 2019;518:18–29. doi:10.1016/j.jcrysgro.2019.04.015. [28]  Tiaden J, Nestler B, Diepers HJ, et al. The multiphase-field model with an integrated concept for modelling solute diffusion. Physica D. 1998;115(1–2):73–86. doi:10.1016/S0167-2789(97)00226-1. [29]  Kim SG, Tae Kim W, Suzuki T, et al. Phase-field modeling of eutectic solidification. J Cryst Growth. 2004;261(1):135–158. doi:10.1016/j.jcrysgro.2003.08.078. [30]  Eiken J, Böttger B, Steinbach I. Multiphase-field approach for multicomponent alloys with extrapolation scheme for numerical application. Phys Rev E. 2006;73(6):066122. doi:10.1103/PhysRevE.73.066122. [31]  Böttger B, Eiken J, Apel M. Multi-ternary extrapolation scheme for efficient coupling of thermodynamic data to a multi-phase-field model. Comput Mater Sci. 2015;108:283–292. doi:10.1016/j.commatsci.2015.03.003. [32]  Morino T, Ode M, Hirosawa S. Direct CALPHAD coupling phase-field model: Closed-form expression for interface composition satisfying equal diffusion potential condition. Phys Rev E. 2024;109(6):065303. doi:10.1103/PhysRevE.109.065303. [33]  Morino T, Ode M, Hirosawa S. An explicit integration approach for predicting the microstructures of multicomponent alloys. Nat Commun. 2025;16(1):6504. doi:10.1038/s41467-025-61246-7.  23  [34]  Kim SG, Kim DI, Kim WT, et al. Computer simulations of two-dimensional and three-dimensional ideal grain growth. Phys Rev E. 2006;74(6):061605. doi:10.1103/PhysRevE.74.061605. [35]  Miyoshi E, Takaki T, Ohno M, et al. Ultra-large-scale phase-field simulation study of ideal grain growth. npj Comput Mater. 2017;3(1):25. doi:10.1038/s41524-017-0029-8. [36]  Miyoshi E, Takaki T, Sakane S, et al. Large-scale phase-field study of anisotropic grain growth: Effects of misorientation-dependent grain boundary energy and mobility. Comput Mater Sci. 2021;186:109992. doi:10.1016/j.commatsci.2020.109992. [37]  Yang Q, Gao Y, Kirshtein A, et al. A free-energy-based and interfacially consistent phase-field model for solid-state sintering without artificial void generation. Comput Mater Sci. 2023;229:112387. doi:10.1016/j.commatsci.2023.112387. [38]  Shi R, Wood M, Heo TW, et al. Towards understanding particle rigid-body motion during solid-state sintering. J Eur Ceram Soc. 2021;41(16):211–231. doi:10.1016/j.jeurceramsoc.2021.09.039. [39]  Ishii A, Kondo K, Yamamoto A, et al. Phase-field modeling of solid-state sintering with interfacial anisotropy. Mater Today Commun. 2023;35:106061. doi:10.1016/j.mtcomm.2023.106061. [40]  Deng J. A Phase Field Model of Sintering with Direction-Dependent Diffusion. Mater Trans. 2012;53(2):385–389. doi:10.2320/matertrans.M2011317. [41]  Biswas S, Schwen D, Wang H, et al. Phase field modeling of sintering: Role of grain orientation and anisotropic properties. Comput Mater Sci. 2018;148:307–319. doi:10.1016/j.commatsci.2018.02.057. [42]  Yang Q, Kirshtein A. A thermodynamically consistent phase-field-micromechanics model of sintering with coupled diffusion and grain motion. J Am Ceram Soc. 2025;108(3):e20279. doi:10.1111/jace.20279. [43]  Chockalingam K, Kouznetsova VG, Van Der Sluis O, et al. 2D Phase field modeling of sintering of silver nanoparticles. Comput Methods Appl Mech Eng. 2016;312:492–508. doi:10.1016/j.cma.2016.07.002. [44]  Andersson J, Ågren J. Models for numerical treatment of multicomponent diffusion in simple phases. J Appl Phys. 1992;72(4):1350–1355. doi:10.1063/1.351745. [45]  Karma A. Phase-Field Formulation for Quantitative Modeling of Alloy Solidification. Phys Rev Lett. 2001;87(11):115701. doi:10.1103/PhysRevLett.87.115701. [46]  Daubner S, Hoffrogge P, Minar M, et al. Triple junction benchmark for multiphase-field and multi-order parameter models. Comput Mater Sci. 2023;219:111995. doi:10.1016/j.commatsci.2022.111995.  25   Figure 1. Profiles of PF variables  and , representing two phases and PEDs. The th material phase is represented by the product .   Figure 2. Calculated phase diagram of a hypothetical A–B–C ternary system at 1273 K, used for model verification in this study.  Table 1 Nondimensionalized parameters used in PF simulations. Parameters Values Surface diffusion mobility,      0.1 Internal interfacial diffusion mobility,    0.01 Volume diffusion mobility,      0.001 Internal interfacial mobility, M 1.0 Spacing of the finite-difference grid,  l 1.0 Thickness of the surface region,  4.0 Width of the internal interfacial region,  4.0 Time increment,  t 0.1   26   Figure 3. Validation of the MPF-PED model in a 2D single-phase system. (a) Initial condition for the simulation: Two circular grains of phase X, possessing the equilibrium concentration, are placed in contact. (b) Resulting state after 2 × 106 computational steps for s / XX = 1.0, showing the formation of a neck. (c) Profiles of  and  ̃  along the white dashed line in (b). The dashed line represents the equilibrium concentration. (d) Comparison of numerical and analytical solutions for dihedral angles. Numerical solutions obtained from simulations with s / XX ranging from 0.2 to 4.0, and analytical solutions calculated using Young’s equation.   Figure 4. Simulation results of 3D multiphase grain growth in a ternary system (phases X, Y, and Z) coexisting with a PED. Snapshots of the morphological evolution for (a)  27  Case 1 (equilibrium initial state) and (b) Case 2 (non-equilibrium initial state for phase X). In (a) and (b), white regions indicate the internal interfaces. (c) Corresponding evolution of the concentration field  ̃A within the material domain for Case 2.   Figure 5. Variation in phase volume fractions (X, Y, and Z) during 3D grain growth simulations. (a) Results for Case 1 (equilibrium initial state). The area enclosed by the dashed square is enlarged. (b) Results for Case 2 (non-equilibrium initial state). Dashed lines represent the equilibrium phase volume fractions calculated using the CALPHAD method based on the overall composition.    28  Statements of Novelty  A novel multi-phase-field model enables thermodynamically consistent simulation of multiphase and multicomponent material processes with vacuum by incorporating a new variable to distinguish material and external domains.   Graphical abstract    29  Supplementary Material  Thermodynamically consistent multi-phase-field model  incorporating passive external domains  Akimitsu Ishiia,*, Yusuke Matsuokaa, Toshiyuki Koyamaa  a Research Center for Structural Materials, National Institute for Materials Science, Tsukuba, Ibaraki 305-0047, Japan *Corresponding author: E-mail address: ISHII.Akimitsu@nims.go.jp  S1. Selection of the interpolation function P() To ensure the model functions correctly, the interpolation function P() is required to continuous and monotonically increasing within the domain 0 ≤  ≤ 1, satisfying boundary conditions P(0) = 0 and P(1) = 1, and possessing point symmetry about  = 0.5. In phase-field modeling, the cubic function P3() = 2(3 − 2) and the quintic P5() = 3(10 − 15 + 62) are commonly employed to meet these criteria. The selection of a specific function is governed by two factors: thermodynamic stability and interface kinetics. First, regarding thermodynamics, the use of a linear function P() =  is insufficient. As shown in Fig. S1, the summation of the double-well potential in Eq. (7) and the double-obstacle potential in Eq. (8) results in a spurious local energy minimum at  ≈ 0.97 at the internal interface ( =  = 0.5) when P() =  is used. Therefore, a polynomial interpolation function with vanishing derivatives at = 0, 1 is necessary to eliminate this artifact and maintain a monotonic energy landscape. Second, regarding kinetics, the choice between P3() and P5() influences the interface mobility. In the time evolution equation for the phase-field variables , P() acts as a coefficient. In the mailto:ISHII.Akimitsu@nims.go.jp 30  region  < 0.5, the derivative of the quintic function is smaller than that of the cubic function (dP()/d < dP()/d). Consequently, the use of P5() results in a stronger suppression of the interface migration velocity near the passive external domain ( ≈ 0), which can lead to numerical stagnation of the interface. Based on these considerations, the cubic function P3() was adopted in the present model.   Figure S1. The barrier energy function fbarrier, which is sum of the double-well potential in Eq. (7) and the double-obstacle potential in Eq. (8), plotted against  at the internal interface condition ( =  = 0.5). The parameters W = 16 and W = 4 were used.  S2. Detailed derivation of Eq. (14) The functional derivative of the total free energy G with respect to  ̃  is given by           chemchemchemi i i iii i i iiP ff f fGc c c ccf f P ff fc c c cc                                               . (S25) Equation (6) yields the following expressions:  31  1ic (S26) and   1ic  . (S27) Therefore, Eq. (S25) can be rewritten as    chemchemi if f P fGf fc c                      . (S28) The last term on the right-hand side of Eq. (S28) can be transformed into chemcNN N ji i i jjf f c fc c c c                  . (S29) Here, ∂   / ∂ ̃  is expressed as follows:  211c cc ccN Nj j jk ki i k k ik kN Nj k j jkikikk k kk kNj jki kkc c cc cc c c c cc c c ccc c cc ccc c                                                  . (S30) where ik is the Kronecker delta. Substituting Eq. (S30) into Eq. (S29) yields the following equation: chemc c cN N NN Nj j jki i j i k jj j kf c f c c fcc c c c c c                               . (S31) According to Eq. (S31), the following expression is derived:  chem chemchem chemc cN NN j jki i i k jj kf f c c ff f cc c c c c                      . (S32) Furthermore, based on the definition of fchem, Eq. (S8) can be rewritten as   32   chemchemc cc cc c cN NN j jki i k jj kN NN N j jki k jj kN N NN j jkj i j kj j kf c c ff cc c c cc c ff cc c cf c f cf cc c c c                                                      , (S33) Substituting Eq. (S33) into Eq. (S28) yields Eq. (14). Alternatively, Eq. (S8) can be transformed as follows:  chemchemchemchemchem chemchemc cc c cccN NN j jki i k jj kN N NN Nj jkj i j kj j kNN Nki kkNki kkf c c ff cc c c cf c f cf cc c c cf ff cc cf ff cc c                                               , (S34) This equation provides the definition of the chemical potential i.    33  S3. Derivation of the partial derivative of    with respect to ci The derivation of the partial derivative ∂   / ∂ci in a multiphase–multicomponent system is presented. The following two vectors are defined: T11 2, , cNc c c  c  (S35) and T11 2, , cNc c c     c . (S36) Using these vectors, the relationship in Eq. (3) can be expressed as follows: N c c . (S37) Partially differentiating both sides of Eq. (S37) with respect to c yields: NN   ccc cI J, (S38) where I is the identity matrix, and J is the Jacobian matrix, defined as 1 11 11 11 1cc ccNN NNc cc cc cc c               J . (S39) Assuming local equilibrium at the interface, the parallel tangent law yields the following relationship:      1 1 2 2 N Nf f f   c c c . (S40) For any two arbitrary phases  and , the following relationship holds:  34           1f ff ff f                                 c cc cc cccc cc c c cH J H JJ H H J, (S41) where H is the Hessian matrix, defined as 2 21 1 112 21 1 11cc c cNN N Nf fc c c cf fc c c c                           H . (S42) Substituting Eq. (S41) into Eq. (S38) yields    1 1N N            I H H J H H J , (S43) Rearranging this equation leads to the following expression: 11 1N        J H H . (S44) The partial derivatives for the first Nc–1 components are determined by calculating the Jacobian matrix via Eq. (S44). The elements for the Ncth component are obtained from the following relationships:  1ccNN ij jic cc c     , (S45) 1ccNi iN jjc ccc    , (S46) and  35  1 1c cccN NN iN ji jc ccc    . (S47)  S4. Preliminary study to determine the outer penalty magnitude To determine the appropriate magnitude for the outer penalty coefficient, simple 2D simulations were performed. As shown in Fig. S2(a), an equilibrium circular phase was placed within the computational domain. The outer potential coefficient was defined as     and simulations were performed for 1 × 105 computational steps under five conditions:    = W, 10W, W, 1000W, and 10000W. Figure S2(b) shows the profile of  along the white dashed line in (a). The case with     = 10000 W resulted in numerical divergence and is therefore excluded from the plot. Conversely, for     = 100W, 1000W, the overshoot of  was sufficiently suppressed. Based on these results, we adopted     = 100W in this study as a value that effectively minimizes overshoot without inducing numerical instability.    Figure S2. Preliminary examination to determine the coefficient of the double-well potential in the regions < 0 and > 1 (a) Initial state for the simulation, (b) Profiles of  36   along the white dashed line in (a) obtained after 1 × 105 computational steps simulation using     = W, 10W, W, and 1000W.