The Allen–Cahn (AC) equation describes how phase separation occurs in binary alloy systems and the dynamics of interfaces between different phases. In the present study, we incorporated the function of high order polynomial potentials in the standard AC equation and present the stability condition for the numerical scheme that is used to solve the AC problem in three-dimensional space. We used a fully explicit Euler scheme for temporal derivatives and a second-order finite difference approach for spatial discretization. However, the explicit scheme is known for its speed and accuracy due to the use of small time steps, but it is subject to a temporal step size limitation. Here, we derived and validated a time step condition that satisfies the discrete maximum principle and assures the stability of the scheme. Several experiments are carried out under the constrained time step to ensure the accuracy of the explicit method, the stability of the scheme, and the discrete maximum principle.

Topics
Linear stability analysis, Phase transitions, Image processing, Alloys, Numerical methods, Computational methods, Computer simulation, Finite difference methods, Mathematical modeling
We study the stability condition for a fully explicit finite difference approach applied to numerically solve the Allen–Cahn (AC) equation in the three-dimensional (3D) domain with high-order polynomial free energy function,15 which is a generalized equation of the conventional AC equation,1 with a zero Neumann boundary condition,
(1)
(2)
where ψ(x, t) represents the order parameter in a domain Ω⊂R3⁠, t represents time, ɛ represents interface width and acquired positive constant value, n is the unit normal vector to the boundary ∂Ω, and Fα(ψ) denotes a double-well free energy function, which is defined as follows:
(3)
Equation (1) can be rewritten by substituting Eq. (3) as follows:
(4)
We note that Eq. (4) can be derived as the L2-gradient flow of the following energy functional:
(5)
Differentiating the energy functional (5) with respect to the time variable and incorporating the zero Neumann boundary condition (2), we can derive the following non-increasing energy property for the governing equation:
ddtE(ψ)=∫ΩψtFα′(ψ)ε2+∇ψt⋅∇ψdx=∫ΩψtFα′(ψ)ε2−Δψdx=−∫Ωψt2dx.

Generally, the AC equation is used for phase separation in materials, solidification processes, and phase transitions between liquid and vapor states. In order to minimize interfacial energy, the evolution of a scalar order parameter is modeled by the AC equation. Various applications of the AC equation include 3D volume reconstruction from slice data, image segmentation,20 shape transformation,5,11 topology optimization,25–27 triply periodic surface reconstruction,17,19 crystal growth,18 and various other problems involving interfacial evolution. The AC equations are typically difficult to solve analytically due to their complexity, so numerical methods are essential for approximating solutions. Lee13 proposed a gradient-decent-like scheme, constructing linearly implicit methods for the AC equation. Ayub et al.2 compared and analyzed operator splitting schemes for solving the AC equation. Zhang et al.29 presented third-order accurate maximum principle preserving schemes for the AC equation. Zhang et al.31 applied high-order strong stability-preserving implicit–explicit Runge–Kutta schemes to numerically solve the space-fractional AC equation, which ensures the preservation of the maximum principle and energy stability. Zhang et al.32 proposed a high-order maximum principle-preserving scheme by employing the Runge–Kutta method for temporal integration and the finite difference method (FDM) for spatial discretization; the study developed a numerical approach. Numerical experiments demonstrated key properties such as the preservation of the maximum principle and energy dissipation. Moreover, a conservative AC equation model and computational methods have been developed. Kwakkel et al.10 redefined the energy functional for the AC equation to prevent mass loss. Xu et al.28 developed modified high-order lattice Boltzmann models for the conservative AC equation. Xia et al.22 presented temporally second-order, unconditionally stable direct schemes for solving the AC and conservative AC equations on the curved surface, which is carried out using a surface mesh made up of piecewise triangular elements and a dual-surface polygonal tessellation.23,24

Higher-order polynomial potential energy allows for the adjustment of the sharpness and propagation speed of phase interfaces,8 enabling the modeling of more complex phase transition dynamics. Compared to traditional fourth-order polynomials, it forms more distinct interfaces and effectively suppresses noise, making it advantageous in applications such as image processing and pattern formation. Furthermore, it can accurately represent intricate interactions, including multi-equilibrium states, making it widely applicable in material science and biological modeling. To solve high-order AC equations, several numerical studies have been conducted. For example, Lee et al.12 employed an interpolation method to address the nonlinear double-well potential term and then applied a second-order operator splitting method to separate the diffusion equation, collectively solving these with Crank–Nicolson and multigrid methods. In addition, Lee et al.16 used effective time step conditions and reformulated the semi-implicit scheme to get a solution of slow dynamics due to high-order polynomial free energy in the AC equation. As we observed in the literature, various numerical methods are applied to solve the high-order AC equation.

Even for unconditionally stable methods, sufficiently small time step sizes are required from the perspective of accuracy. In this context, while the fully explicit scheme has a time step constraint, it remains one of the fastest and most accurate methods for solving the AC equation. In addition, it can effectively be applied in boundary problems in irregular domains and adaptive domains. In particular, high-order potential enables the investigation of complex 3D interfacial problems. Therefore, stability analysis of the explicit method for the high-order AC equation is a crucial topic of study. More recently, Choi et al.3 employed a fully explicit finite difference scheme to solve a two-dimensional (2D) high-order AC model, successfully ensuring the preservation of the maximum principle. However, examination of the stability condition of the 3D AC model using a fully explicit scheme with a high-order polynomial potential has not yet been explored. Consequently, this paper aims to investigate the stability conditions of the 3D high-order AC equation.

The section-wise organization of the paper is as follows: Sec. II presents a fully explicit numerical solution and analyzes the stability of the high-order AC equation. A numerical proof of the proposed theorem is provided in Sec. III. As a final note, Sec. IV concludes.

In this study, we propose a temporal step restriction to guarantee the maximum principle’s compliance. In addition, we used a fully explicit method for the high-order AC equation. Let the AC equation be defined over a 3D domain Ω = (Xl, Xr) × (Yl, Yr) × (Zl, Zr) with the grid size h = (Xr − Xl)/Nx = (Yr − Yl)/Ny = (Zr − Zl)/Nz. Here, h is the spatial step size, and Nx, Ny, and Nz are the number of spatial grids in the x-, y-, and z-directions, respectively. In addition, the discrete domain is defined as Ωh = {(xa, yb, zc)|xa = Xl + (a − 0.5)h, yb = Yl + (b − 0.5)h, zc = Zl + (c − 0.5)h a = 1, 2, …, Nx, b = 1, 2, …, Ny, c = 1, 2, …, Nz}. On the computational domain, ψ(xa, yb, zc, nΔt) is approximated as ψabcn for n = 1, 2, …, Nt, where Δt = T/Nt with final time T and number of time steps Nt. We define the discrete maximum norm by ‖ψn‖∞=max1≤i≤Nx,1≤j≤Ny,1≤k≤Nz|ψijkn|⁠,
(6)
We may use a 27-point stencil isotropic finite difference scheme14 that can be utilized for isotropic solutions. In addition, homogeneous Neumann boundary conditions are applied within the 3D computational domain for a = 1, 2, …, Nx, b = 1, 2, …, Ny, c = 1, 2, …, Nz, as follows:
ψ0,bc=ψ1,bc,ψNx+1,bc=ψNx,bc,ψa,0,c=ψa,1,c,ψa,Ny+1,c=ψa,Ny,c,ψab,0=ψab,1,ψab,Nz+1=ψab,Nz.

Next, for the fully explicit FDM approach for the 3D high-order AC equation, we shall prove a theorem to meet the maximum principle.

Theorem 1.
Under the assumption that ‖ψ0‖∞ ≤ 1 is satisfied by the initial condition, the fully explicit numerical technique outlined in (6) guarantees that the numerical solutions stay bounded at every discrete time step,
‖ψn+1‖∞≤1,forn≥0,
if
(7)

Proof.
It is assumed that the maximum norm of the initial state, denoted as ‖ψ0‖∞, is less than or equal to 1. Using mathematical induction, we assume that the maximum norm of the state at the nth step, denoted as ‖ψn‖∞, is also less than or equal to 1. We then express the explicit computational method for the AC model in the 3D domain, as described by (6), in the following manner:
(8)
From (8), since ‖ψn‖∞ ≤ 1, the following inequalities can be derived:
(9)
and
(10)

Our goal is to identify the criteria for Δt to ensure that LB ≥ −1 and UB ≤ 1, which concludes ‖ψn+1‖∞ ≤ 1.

Case 1.

LB ≥ −1

From the inequality (9), it can be reformulated as
(11)
It can be reformulated as
(12)
If ψabcn=−1⁠, the time step is unconstrained. Otherwise, when ψabcn>−1⁠, dividing both sides of (12) by the negative quantity −(1+ψjin)<0⁠, we obtain
(13)
If −1≤ψabcn≤1⁠, then given that 0≤1−(ψabcn)α1−ψabcn+⋯+(−1)α−1(ψabcn)α−1≤2α⁠, we have
(14)
As a result, we have the following condition:
(15)

Next, it is necessary to determine the condition of Δt that ensures

Case 2.

UB ≤ 1

From the inequality (10), the following inequality can be restated:
(16)
If ψabcn=1⁠, then Δt is unconstrained. When ψabcn<1⁠, dividing both sides of (16) by (1−ψabcn)>0 yields
Δt(ψabcn)2α−11+ψabcn+⋯+(ψabcn)2α−1αε2+6h2≤1.
Therefore, Δt must satisfy the following restriction:
Δt≤αε2h2h2(ψabcn)2α−11+ψabcn+⋯+(ψabcn)2α−1+6αε2.
Given ψabcn≤1⁠, we can obtain the following condition:
Δt≤ε2h22h2+3ε2.

Hence, Δt≤0.5ε2h2/h2+2ε2 ensure that ψabcn+1≤1⁠.□

Remark 1.

In the field of numerical analysis, demonstrating the stability of numerical schemes is often a challenging task, especially when directly proving the boundedness of numerical solutions. In such cases, energy stability has been commonly used as an alternative to establish the stability of the numerical method.7, However, the proposed time step constraint in this study ensures the maximum bounded principle for the numerical solution in the explicit method, thereby directly proving the stability of the numerical solution without the need for additional energy stability analysis. Nevertheless, it is worth noting that if the proposed time step constraint is satisfied, energy stability may still be demonstrated by applying the approach outlined in the references for one-dimensional space4  and higher-dimensional spaces.21 

We note that4 it established the situation when α = 1. In this work, we use mathematical induction to extend the stability analysis of numerical approaches for the high-order AC model in a 3D domain. Furthermore, Lemma 3.7 from Ref. 30 may be used to evaluate this time step limitation. Subsequently, we prove that the temporal step size requirement specified in Theorem 1 is the best feasible, since it reflects the largest temporal step that follows the maximal principle. Let
Δtmax=ε2h22(h2+3ε2),
be the maximum number of time steps required to fulfill the maximum principle for the high-order AC equation. Assume Δt > Δtmax in order to verify that Δtmax is optimum. This implies 0 < Δtmax/Δt < 1. Examine the following polynomial function that is defined on the (0, 1) interval:
f(ϕ)=6αε2ϕ+h2ϕ2α(1+ϕ+ϕ2+⋯+ϕ2α−1)2α(h2+2ε2).
Given that f′(ψ) > 0 for ψ ∈ (0, 1), consequently the function above must increase monotonically. Given that f(ψ) is continuous, f(0) = 0, and f(1) = 1, there exists a unique value of ψ such that Δtmax/Δt = f(ψ). Let
ψ(xa,yb,zc,0)=ϕ, ifa=p,b=q,c=r,1, otherwise,
where 1 < p < Nx, 1 < q < Ny, and 1 < r < Nz. From Eq. (8), we get
ψpqr1=ψpqr0+Δt(ψpqr0)2α−11−(ψpqr0)2ααε2+ψp+1,q,r0+ψp−1,q,r0+ψp,q+1,r0+ψp,q−1,r0+ψp,q,r−10+ψp,q,r+10−6ψpqr0h2=ϕ+Δtmaxf(ϕ)ϕ2α−11−ϕ2ααε2+6−6ϕh2=ϕ+αε2h26αε2ϕ+h2ϕ2α(1+ϕ+⋯+ϕ2α−1)ϕ2α−11−ϕ2ααε2+6−6ϕh2=ϕ+1ϕ−1>1.
This implies that for any temporal step sizes greater than Δtmax, the maximum principle will not be maintained. Δtmax is the most suitable bound as a consequence. We note that the maximum principle is important, and some numerical methods do not satisfy this property.9 For future work, we will apply the current analysis to the AC equation on a triangular mesh.6 

We conduct numerical experiments to verify the maximum principle, analyze mean curvature flow, and examine how the parameter α impacts the equation’s dynamic behavior.

Prior to performing numerical simulations, we evaluate the accuracy of the numerical methodology. For this purpose, we use the initial condition given below:
(17)
The parameters employed to evaluate the spatial accuracy are Nx = 2n for n = 5, 6, 7, 8, and Nz = Ny = Nx on Ω = (−1, 1) × (−1, 1) × (−1, 1), ɛ = 0.001, Δt = 10−7, for a maximum of 100 iterations, while for the reference value we used Nx = 24 and h = 1/8. Then, we calculated the discrete l2-norm error in terms of relative error for different values of α = 4, 6, 8, 10, 12. The discrete l2-norm ‖·‖2 is defined as follows:
‖ψ‖2=1NxNyNz∑k=1Nz∑j=1Ny∑i=1Nx(ψijk−ψref)2,
where ψijk is the value at spatial step size h, and ψref is the reference value at spatial step size h/2. Subsequently, we determined the spatial convergence rate, defined as log2(‖E1‖2/‖E2‖2), where E1 denotes the error with a spatial step size h and E2 denotes the error with a spatial step size h/2. The spatial error with convergence rate is shown in Table I.
TABLE I.

The space-dependent error and its rate of convergence.

h = 1/32h = 1/64h = 1/128h = 1/256
α = 4 Error rate 9.9495 × 10−3 1.958 2.5614 × 10−3 1.994 0.6432 × 10−3 1.998 0.1610 × 10−3 
α = 6 Error rate 9.9849 × 10−3 1.974 2.5408 × 10−3 1.994 0.6379 × 10−3 1.998 0.1596 × 10−3 
α = 8 Error rate 9.9914 × 10−3 1.976 2.5404 × 10−3 1.994 0.6378 × 10−3 1.998 0.1596 × 10−3 
α = 10 Error rate 9.9918 × 10−3 1.976 2.5404 × 10−3 1.994 0.6378 × 10−3 1.998 0.1596 × 10−3 
α = 12 Error rate 9.9918 × 10−3 1.976 2.5404 × 10−3 1.994 0.6378 × 10−3 1.998 0.1596 × 10−3 
h = 1/32h = 1/64h = 1/128h = 1/256
α = 4 Error rate 9.9495 × 10−3 1.958 2.5614 × 10−3 1.994 0.6432 × 10−3 1.998 0.1610 × 10−3 
α = 6 Error rate 9.9849 × 10−3 1.974 2.5408 × 10−3 1.994 0.6379 × 10−3 1.998 0.1596 × 10−3 
α = 8 Error rate 9.9914 × 10−3 1.976 2.5404 × 10−3 1.994 0.6378 × 10−3 1.998 0.1596 × 10−3 
α = 10 Error rate 9.9918 × 10−3 1.976 2.5404 × 10−3 1.994 0.6378 × 10−3 1.998 0.1596 × 10−3 
α = 12 Error rate 9.9918 × 10−3 1.976 2.5404 × 10−3 1.994 0.6378 × 10−3 1.998 0.1596 × 10−3 
View Large

Subsequently, we evaluated the temporal accuracy of the method for Δtmax = 10−4 using the same initial condition as specified in (17), with a perturbation value of 0.4 and the following parameters: Nx = Ny = Nz = 27, Δt = Δtmax/2n for n = 0, 1, 2, 3, 4, ɛ = h/10, and T = 100Δtmax, including α = 4, 6, 8, 10, 12. Here for reference value we used Δt = Δtmax. Table II presents time-dependent error and the corresponding convergence rate. In addition, the graphical representation of spatial and temporal error is present in Fig. 1.

TABLE II.

Time-dependent error and its rate of convergence.

Δt = Δtmax/2Δt = Δtmax/4Δt = Δtmax/8Δt = Δtmax/16
α = 4 Error rate 7.3140 × 10−4 1 3.6570 × 10−4 1 1.8285 × 10−4 1 0.9143 × 10−4 
α = 6 Error rate 8.3219 × 10−4 1 4.1609 × 10−4 1 2.0804 × 10−4 1 1.0402 × 10−4 
α = 8 Error rate 8.3673 × 10−4 1 4.1837 × 10−4 1 2.0918 × 10−4 1 1.0459 × 10−4 
α = 10 Error rate 8.3679 × 10−4 1 4.1840 × 10−4 1 2.0920 × 10−4 1 1.0460 × 10−4 
α = 12 Error rate 8.3679 × 10−4 1 4.1840 × 10−4 1 2.0920 × 10−4 1 1.0460 × 10−4 
Δt = Δtmax/2Δt = Δtmax/4Δt = Δtmax/8Δt = Δtmax/16
α = 4 Error rate 7.3140 × 10−4 1 3.6570 × 10−4 1 1.8285 × 10−4 1 0.9143 × 10−4 
α = 6 Error rate 8.3219 × 10−4 1 4.1609 × 10−4 1 2.0804 × 10−4 1 1.0402 × 10−4 
α = 8 Error rate 8.3673 × 10−4 1 4.1837 × 10−4 1 2.0918 × 10−4 1 1.0459 × 10−4 
α = 10 Error rate 8.3679 × 10−4 1 4.1840 × 10−4 1 2.0920 × 10−4 1 1.0460 × 10−4 
α = 12 Error rate 8.3679 × 10−4 1 4.1840 × 10−4 1 2.0920 × 10−4 1 1.0460 × 10−4 
View Large
FIG. 1.
Graphical representation of (a) spatial and (b) temporal error.
View largeDownload slide

Graphical representation of (a) spatial and (b) temporal error.

FIG. 1.
Graphical representation of (a) spatial and (b) temporal error.
View largeDownload slide

Graphical representation of (a) spatial and (b) temporal error.

Close modal

Therefore, it was concluded that the scheme demonstrates a second-order spatial error rate and a first-order temporal error rate.

Now, we check that our algorithm preserves the maximum principle with t = Δtmax. For this purpose, we used the following initial condition:
ψ(x,y,0)=0.7⁡cos(2πx)cos(2πy)cos(2πz),
on Ω = (0, 1) × (0, 1) × (0, 1) with Nx = 50 = Ny = Nz, α = 3, h = 0.02, ɛ = h, Δtmax = ɛ2h2/(2(h2 + 3ɛ2)), and the final time T = 0.06 with maximum iterations = T/Δtmax. For Δt = 1.18tmax, we calculated whether our numerical scheme preserves the maximum principle or not and found that our scheme preserves the maximum principle. Figure 2 shows the time evolution of the maximum and minimum values of ψ for the time steps Δt = Δtmax and Δt = 1.18Δtmax. If a time step larger than Δtmax is used, the value of ψ is observed to go beyond the limits of −1 and 1.
FIG. 2.
Temporal evolutions of the maximum and minimum values of ψn for the time steps Δt = Δtmax and Δt = 1.18Δtmax.
View largeDownload slide

Temporal evolutions of the maximum and minimum values of ψn for the time steps Δt = Δtmax and Δt = 1.18Δtmax.

FIG. 2.
Temporal evolutions of the maximum and minimum values of ψn for the time steps Δt = Δtmax and Δt = 1.18Δtmax.
View largeDownload slide

Temporal evolutions of the maximum and minimum values of ψn for the time steps Δt = Δtmax and Δt = 1.18Δtmax.

Close modal
In this subsection, we investigate the energy variation of the computational solution to the governing equation. For the computational solution of the AC equation using a fully explicit method, it is shown that the discrete energy dissipation law is satisfied, provided the time step condition derived from the stability analysis is met. Here, the discrete energy functional with the homogeneous Neumann boundary condition is defined as follows:
E3(ψn)=h3∑i=1Nx∑j=1Ny∑k=1NzFα(ψijkn)ε2+h2∑j=1Ny∑k=1Nz∑i=1Nx−1(ψi+1,j,kn−ψi,j,kn)2+h2∑i=1Nx∑k=1Nz∑j=1Ny−1(ψi,j+1,kn−ψi,j,kn)2+h2∑i=1Nx∑j=1Ny∑j=1Nz−1(ψi,j,k+1n−ψi,j,kn)2.
Here, α = 3, Δt = Δtmax is chosen, h = 1/64 is uniform mesh size, and ɛ = h in the 3D unit domain (0, 1) × (0, 1) × (0, 1). The random perturbed initial condition is given as follows:
ψ(x,y,z,0)=0.1×rand(x,y,z),(x,y,z)∈(0,1)×(0,1)×(0,1),
where rand(x, y, z) is a uniformly distributed random number between −1 and 1. Figure 3 shows the time evolution of the zero-level isosurface of the numerical solutions and the discrete energy functional. It demonstrates that the proposed time step constraint ensures the energy dissipation law of the high-order AC equation.
FIG. 3.
Snapshots of the zero-level isosurface of the computational solutions with temporal evolution of discrete energy functional.
View largeDownload slide

Snapshots of the zero-level isosurface of the computational solutions with temporal evolution of discrete energy functional.

FIG. 3.
Snapshots of the zero-level isosurface of the computational solutions with temporal evolution of discrete energy functional.
View largeDownload slide

Snapshots of the zero-level isosurface of the computational solutions with temporal evolution of discrete energy functional.

Close modal
In the AC equation, motion by mean curvature describes how phase boundaries evolve over time. This deals with the dynamics of the interface adjusting to reduce the total interface energy between two different phases. The interface movement is governed by the mean curvature, which works to eliminate imperfections and decrease surface area, controlling how the interface moves. Numerous physical processes, such as phase separation, pattern creation, and the generation of material microstructures, depend on this phenomenon. In this work, we use a completely explicit method for the 3D AC equation with a high-order polynomial potential in order to analyze motion by mean curvature. Let
ψ(x,y,z,0)=tanhR0−x2+y2+z22ε.
We used the following parameters: Nx = Ny = Nz = 128, α = 3, R0 = 0.8, and ε=mh/22tanh(0.9)⁠, where m = 3.8, Δt = h2ɛ2/(2h2 + 6ɛ2), and simulation time T = 0.15 for which Nt = (T/Δt). The exact radius is defined as R(t)=r2−4t in 3D space. The time evolution of the zero-level isosurface of the 3D high-order AC equation is shown in Fig. 4, along with a comparison of the analytical and numerical findings. Together, the results from Fig. 4 show that, in the case of the explicit method, the studied time step guarantees method stability, and the solutions correctly follow the mean curvature flow. The computational experiment demonstrates a high degree of agreement between the numerical findings and the theoretical radius, with the numerical results aligning closely with the two.
FIG. 4.
Comparison between analytical and numerical radius with the temporal progression of the computational results illustrated as an isosurface when Δt = h2ɛ2/(2h2 + 6ɛ2).
View largeDownload slide

Comparison between analytical and numerical radius with the temporal progression of the computational results illustrated as an isosurface when Δt = h2ɛ2/(2h2 + 6ɛ2).

FIG. 4.
Comparison between analytical and numerical radius with the temporal progression of the computational results illustrated as an isosurface when Δt = h2ɛ2/(2h2 + 6ɛ2).
View largeDownload slide

Comparison between analytical and numerical radius with the temporal progression of the computational results illustrated as an isosurface when Δt = h2ɛ2/(2h2 + 6ɛ2).

Close modal

We analyzed the stability and accuracy of an explicit computational method for solving the high-order AC equation in three-dimensional space, essential for modeling complex phase dynamics in alloy systems. By incorporating high-order polynomial potentials, we addressed the requirements for capturing stable and metastable phases, anisotropy, and intricate interface behaviors. Using an explicit Euler scheme for time discretization and a second-order FDM for spatial discretization, we derived a critical time step condition that ensures the discrete maximum principle and stabilizes the scheme. Our experiments validated that, under this constrained time step, the explicit scheme maintains stability, accuracy, andadherence to the discrete maximum principle. This work provides a reliable approach to capturing complex phase dynamics in the AC models, contributing to improved numerical simulations in multi-phase material systems.

The corresponding author (J. S Kim) was supported by Korea University Grant. Jyoti was supported by the Brain Pool program funded by the Ministry of Science and ICT through the National Research Foundation of Korea (Grant No. 2022H1D3A2A02081237). We appreciate the reviewers for their insightful comments and constructive suggestions, which have greatly improved the quality of our paper.

The authors have no conflicts to disclose.

Seokjun Ham: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Jyoti: Formal analysis (equal); Investigation (equal); Methodology (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Jaeyong Choi: Formal analysis (equal); Investigation (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Yunjae Nam: Software (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Junseok Kim: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article.

The MATLAB code is provided as follows:

% Motion by mean curvature in 3D space 
clear; clear all; close all; clc; 
Nx = 128; xl = -1; xr = 1; yl = -1; yr = 1; zl = -1; zr = 1; h = (xr-xl)/Nx; 
x = linspace(xl-0.5*h,xr+0.5*h,Nx+2); y = linspace(yl-0.5*h,yr+0.5*h,Nx+2); 
z = linspace(zl-0.5*h,zr+0.5*h,Nx+2);[xx, yy, zz] = ndgrid(x, y, z); 
psi =zeros(Nx+2,Nx+2,Nx+2); 
alpha = 3; m=3.8; ep=(m*h)/(2*sqrt(2)*atanh(0.9)); ep2=alpha*ep^2; 
dt= (h^2*ep^2)/(2*h^2+6*ep^2); 
for i=1:Nx+2 
for j=1:Nx+2 
for k=1:Nx+2 
psi(i,j,k) = (tanh((0.7-sqrt((x(i)).^2+(y(j)).^2+(z(k))^2))/(sqrt(2)*ep))); 
end 
end 
end 
figure (1); clf; grid on; box on; hold on; 
set(gca,′FontSize′,18) 
fv=isosurface(x,y,z,permute(psi,[2 1 3]),0); c=fv.vertices; is=patch(fv); 
camlight headlight; is.FaceColor=′[0 1 1]′; is.FaceAlpha=1; 
is.EdgeColor=′none′; 
axis image; box on; view(3); axis([xl xr yl yr zl zr]) 
text(′String′,′$x$′,′Position′,[0.877 -0.778 -1.436], ... 
′interpreter′,′latex′,′fontsize′,20) 
text(′String′,′$y$′,′Position′,[-1.057 0.842 -1.249], ... 
′interpreter′,′latex′,′fontsize′,20) 
text(′String′,′$z$′,′Position′,[-0.688 1.711 0.326], ... 
′interpreter′,′latex′,′fontsize′,20) 
print(′3D_sphere_initial.eps′,′-depsc′) 
radi(1)= mean(sqrt(c(:,1).^2+c(:,2).^2+c(:,3).^2)); 
ii=2:Nx+1; jj=2:Nx+1; kk=2:Nx+1; T = 0.1; Nt=round(T/dt); 
for iter = 1:Nt 
psi(1,:,:)=psi(2,:,:); psi(end,:,:)=psi(end-1,:,:); 
psi(:,1,:)=psi(:,2,:); psi(:,end,:)=psi(:,end-1,:); 
psi(:,:,1)=psi(:,:,2); psi(:,:,end)=psi(:,:,end-1); 
psi(ii,jj,kk) = psi(ii,jj,kk) ... 
+dt*(psi(ii,jj,kk).^(2*alpha-1).*(1-psi(ii,jj,kk).^(2*alpha))/ep2 ... 
+(psi(ii+1,jj,kk)+psi(ii-1,jj,kk)+psi(ii,jj+1,kk)+psi(ii,jj-1,kk) ... 
+psi(ii,jj,kk+1)+psi(ii,jj,kk-1)-6*psi(ii,jj,kk))/h^2); 
fv=isosurface(x,y,z,psi,0); c=fv.vertices; 
radi(iter+1)= mean(sqrt(c(:,1).^2+c(:,2).^2+c(:,3).^2)); 
if abs(iter*dt-0.08)<0.7*dt 
figure (2); clf; hold on; box on; grid on; 
set(gca,′FontSize′,18) 
fv=isosurface(x,y,z,permute(psi,[2 1 3]),0); c=fv.vertices; 
is=patch(fv); camlight headlight; is.FaceColor=′[0 1 1]′; is.FaceAlpha=1; 
is.EdgeColor=′none′;axis image; box on; view(3); axis([xl xr yl yr zl zr]) 
text(′String′,′$x$′,′Position′,[0.877 -0.778 -1.436], ... 
′interpreter′,′latex′,′fontsize′,20) 
text(′String′,′$y$′,′Position′,[-1.057 0.842 -1.249], ... 
′interpreter′,′latex′,′fontsize′,20) 
text(′String′,′$z$′,′Position′,[-0.688 1.711 0.326], ... 
′interpreter′,′latex′,′fontsize′,20) 
print(′3D_sphere_initial.eps′,′-depsc′) 
pause(0.05) 
print(′3D_sphere_1.eps′,′-depsc′) 
end 
end 
figure (3); clf; hold on; box on; grid on; 
set(gca,′FontSize′,18);fv=isosurface(x,y,z,permute(psi,[2 1 3]),0);c=fv.vertices; 
is=patch(fv); camlight headlight; is.FaceColor=′[0 1 1]′; is.FaceAlpha=1; 
is.EdgeColor=′none′;axis image; box on; view(3); axis([xl xr yl yr zl zr]); 
text(′String′,′$x$′,′Position′,[0.877 -0.778 -1.436], ... 
′interpreter′,′latex′,′fontsize′,20) 
text(′String′,′$y$′,′Position′,[-1.057 0.842 -1.249], ... 
′interpreter′,′latex′,′fontsize′,20) 
text(′String′,′$z$′,′Position′,[-0.688 1.711 0.326], ... 
′interpreter′,′latex′,′fontsize′,20) 
print(′3D_sphere_initial.eps′,′-depsc′) 
pause(0.05);print(′3D_sphere_2.eps′,′-depsc′) 
figure (4); clf; hold on; box on; 
set(gca,′FontSize′,18); set(gcf,′Position′,[397.3680981595091 548.7300613496932 ... 
770.9447852760734 342.1840490797546]); tt = linspace(0,Nt*dt,Nt+1); 
plot(tt,sqrt(0.7^2-4*(tt)),′k-′,′LineWidth′,1.2) 
plot(tt,radi,′ko′,′MarkerIndices′,1:185:Nt,′MarkerSize′,8,′LineWidth′,1.2) 
axis([0 0.1 0.3 0.7]);text(′String′,′$t$′,′Position′,[0.0929 0.264], ... 
′interpreter′,′latex′,′fontsize′,20); 
text(′String′,′$R(t)$′,′Position′,[-0.009 0.653], ... 
′interpreter′,′latex′,′fontsize′,20); 
legend(′Analytic result′,′Numerical result′);print(′meancurv_3d.eps′,′-depsc′) 
% Motion by mean curvature in 3D space 
clear; clear all; close all; clc; 
Nx = 128; xl = -1; xr = 1; yl = -1; yr = 1; zl = -1; zr = 1; h = (xr-xl)/Nx; 
x = linspace(xl-0.5*h,xr+0.5*h,Nx+2); y = linspace(yl-0.5*h,yr+0.5*h,Nx+2); 
z = linspace(zl-0.5*h,zr+0.5*h,Nx+2);[xx, yy, zz] = ndgrid(x, y, z); 
psi =zeros(Nx+2,Nx+2,Nx+2); 
alpha = 3; m=3.8; ep=(m*h)/(2*sqrt(2)*atanh(0.9)); ep2=alpha*ep^2; 
dt= (h^2*ep^2)/(2*h^2+6*ep^2); 
for i=1:Nx+2 
for j=1:Nx+2 
for k=1:Nx+2 
psi(i,j,k) = (tanh((0.7-sqrt((x(i)).^2+(y(j)).^2+(z(k))^2))/(sqrt(2)*ep))); 
end 
end 
end 
figure (1); clf; grid on; box on; hold on; 
set(gca,′FontSize′,18) 
fv=isosurface(x,y,z,permute(psi,[2 1 3]),0); c=fv.vertices; is=patch(fv); 
camlight headlight; is.FaceColor=′[0 1 1]′; is.FaceAlpha=1; 
is.EdgeColor=′none′; 
axis image; box on; view(3); axis([xl xr yl yr zl zr]) 
text(′String′,′$x$′,′Position′,[0.877 -0.778 -1.436], ... 
′interpreter′,′latex′,′fontsize′,20) 
text(′String′,′$y$′,′Position′,[-1.057 0.842 -1.249], ... 
′interpreter′,′latex′,′fontsize′,20) 
text(′String′,′$z$′,′Position′,[-0.688 1.711 0.326], ... 
′interpreter′,′latex′,′fontsize′,20) 
print(′3D_sphere_initial.eps′,′-depsc′) 
radi(1)= mean(sqrt(c(:,1).^2+c(:,2).^2+c(:,3).^2)); 
ii=2:Nx+1; jj=2:Nx+1; kk=2:Nx+1; T = 0.1; Nt=round(T/dt); 
for iter = 1:Nt 
psi(1,:,:)=psi(2,:,:); psi(end,:,:)=psi(end-1,:,:); 
psi(:,1,:)=psi(:,2,:); psi(:,end,:)=psi(:,end-1,:); 
psi(:,:,1)=psi(:,:,2); psi(:,:,end)=psi(:,:,end-1); 
psi(ii,jj,kk) = psi(ii,jj,kk) ... 
+dt*(psi(ii,jj,kk).^(2*alpha-1).*(1-psi(ii,jj,kk).^(2*alpha))/ep2 ... 
+(psi(ii+1,jj,kk)+psi(ii-1,jj,kk)+psi(ii,jj+1,kk)+psi(ii,jj-1,kk) ... 
+psi(ii,jj,kk+1)+psi(ii,jj,kk-1)-6*psi(ii,jj,kk))/h^2); 
fv=isosurface(x,y,z,psi,0); c=fv.vertices; 
radi(iter+1)= mean(sqrt(c(:,1).^2+c(:,2).^2+c(:,3).^2)); 
if abs(iter*dt-0.08)<0.7*dt 
figure (2); clf; hold on; box on; grid on; 
set(gca,′FontSize′,18) 
fv=isosurface(x,y,z,permute(psi,[2 1 3]),0); c=fv.vertices; 
is=patch(fv); camlight headlight; is.FaceColor=′[0 1 1]′; is.FaceAlpha=1; 
is.EdgeColor=′none′;axis image; box on; view(3); axis([xl xr yl yr zl zr]) 
text(′String′,′$x$′,′Position′,[0.877 -0.778 -1.436], ... 
′interpreter′,′latex′,′fontsize′,20) 
text(′String′,′$y$′,′Position′,[-1.057 0.842 -1.249], ... 
′interpreter′,′latex′,′fontsize′,20) 
text(′String′,′$z$′,′Position′,[-0.688 1.711 0.326], ... 
′interpreter′,′latex′,′fontsize′,20) 
print(′3D_sphere_initial.eps′,′-depsc′) 
pause(0.05) 
print(′3D_sphere_1.eps′,′-depsc′) 
end 
end 
figure (3); clf; hold on; box on; grid on; 
set(gca,′FontSize′,18);fv=isosurface(x,y,z,permute(psi,[2 1 3]),0);c=fv.vertices; 
is=patch(fv); camlight headlight; is.FaceColor=′[0 1 1]′; is.FaceAlpha=1; 
is.EdgeColor=′none′;axis image; box on; view(3); axis([xl xr yl yr zl zr]); 
text(′String′,′$x$′,′Position′,[0.877 -0.778 -1.436], ... 
′interpreter′,′latex′,′fontsize′,20) 
text(′String′,′$y$′,′Position′,[-1.057 0.842 -1.249], ... 
′interpreter′,′latex′,′fontsize′,20) 
text(′String′,′$z$′,′Position′,[-0.688 1.711 0.326], ... 
′interpreter′,′latex′,′fontsize′,20) 
print(′3D_sphere_initial.eps′,′-depsc′) 
pause(0.05);print(′3D_sphere_2.eps′,′-depsc′) 
figure (4); clf; hold on; box on; 
set(gca,′FontSize′,18); set(gcf,′Position′,[397.3680981595091 548.7300613496932 ... 
770.9447852760734 342.1840490797546]); tt = linspace(0,Nt*dt,Nt+1); 
plot(tt,sqrt(0.7^2-4*(tt)),′k-′,′LineWidth′,1.2) 
plot(tt,radi,′ko′,′MarkerIndices′,1:185:Nt,′MarkerSize′,8,′LineWidth′,1.2) 
axis([0 0.1 0.3 0.7]);text(′String′,′$t$′,′Position′,[0.0929 0.264], ... 
′interpreter′,′latex′,′fontsize′,20); 
text(′String′,′$R(t)$′,′Position′,[-0.009 0.653], ... 
′interpreter′,′latex′,′fontsize′,20); 
legend(′Analytic result′,′Numerical result′);print(′meancurv_3d.eps′,′-depsc′) 
View Large

1.
S. M.
Allen
 and
J. W.
Cahn
, “
A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening
,”
Acta Metall.
27
,
1085
–
1095
(
1979
).
https://doi.org/10.1016/0001-6160(79)90196-2
Google Scholar
Crossref
Search ADS
 
2.
S.
Ayub
,
H.
Affan
, and
A.
Shah
, “
Comparison of operator splitting schemes for the numerical solution of the Allen–Cahn equation
,”
AIP Adv.
9
(
12
),
125202
(
2019
).
https://doi.org/10.1063/1.5126651
Google Scholar
Crossref
Search ADS
 
3.
J.
Choi
,
S.
Ham
,
S.
Kwak
,
Y.
Hwang
, and
J.
Kim
, “
Stability analysis of an explicit numerical scheme for the Allen–Cahn equation with high-order polynomial potentials
,”
AIMS Math.
9
(
7
),
19332
–
19344
(
2024
).
https://doi.org/10.3934/math.2024941
Google Scholar
Crossref
Search ADS
 
4.
S.
Ham
 and
J.
Kim
, “
Stability analysis for a maximum principle preserving explicit scheme of the Allen–Cahn equation
,”
Math. Comput. Simul.
207
,
453
–
465
(
2023
).
https://doi.org/10.1016/j.matcom.2023.01.016
Google Scholar
Crossref
Search ADS
 
5.
Z.
Han
,
H.
Xu
, and
J.
Wang
, “
A simple shape transformation method based on phase-field model
,”
Comput. Math. Appl.
147
,
121
–
129
(
2023
).
https://doi.org/10.1016/j.camwa.2023.07.020
Google Scholar
Crossref
Search ADS
 
6.
Y.
Hwang
,
S.
Ham
,
C.
Lee
,
G.
Lee
,
S.
Kang
, and
J.
Kim
, “
A simple and efficient numerical method for the Allen–Cahn equation on effective symmetric triangular meshes
,”
Electron. Res. Arch.
31
(
8
),
4557
–
4578
(
2023
).
https://doi.org/10.3934/era.2023233
Google Scholar
Crossref
Search ADS
 
7.
J.
Kim
, “
Phase-field models for multi-component fluid flows
,”
Commun. Comput. Phys.
12
(
3
),
613
–
661
(
2012
).
https://doi.org/10.4208/cicp.301110.040811a
Google Scholar
Crossref
Search ADS
 
8.
J.
Kim
, “
Modified wave-front propagation and dynamics coming from higher-order double-well potentials in the Allen–Cahn equations
,”
Mathematics
12
(
23
),
3796
(
2024
).
https://doi.org/10.3390/math12233796
Google Scholar
Crossref
Search ADS
 
9.
J.
Kim
,
S.
Kwak
,
H. G.
Lee
,
Y.
Hwang
, and
S.
Ham
, “
A maximum principle of the Fourier spectral method for diffusion equations
,”
Electron. Res. Arch.
31
(
9
),
5396
–
5405
(
2023
).
https://doi.org/10.3934/era.2023273
Google Scholar
Crossref
Search ADS
 
10.
M.
Kwakkel
,
M.
Fernandino
, and
C. A.
Dorao
, “
A redefined energy functional to prevent mass loss in phase-field methods
,”
AIP Adv.
10
(
6
),
065124
(
2020
).
https://doi.org/10.1063/1.5142353
Google Scholar
Crossref
Search ADS
 
11.
S.
Lai
,
B.
Jiang
,
Q.
Xia
,
B.
Xia
,
J.
Kim
, and
Y.
Li
, “
On the phase-field algorithm for distinguishing connected regions in digital model
,”
Eng. Anal. Boundary Elem.
168
,
105918
(
2024
).
https://doi.org/10.1016/j.enganabound.2024.105918
Google Scholar
Crossref
Search ADS
 
12.
C.
Lee
,
H.
Kim
,
S.
Yoon
,
S.
Kim
,
D.
Lee
,
J.
Park
, and
J.
Kim
, “
An unconditionally stable scheme for the Allen–Cahn equation with high-order polynomial free energy
,”
Commun. Nonlinear Sci. Numer. Simul.
95
,
105658
(
2021
).
https://doi.org/10.1016/j.cnsns.2020.105658
Google Scholar
Crossref
Search ADS
 
13.
D.
Lee
, “
Gradient-descent-like scheme for the Allen–Cahn equation
,”
AIP Adv.
13
(
8
),
085010
(
2023
).
https://doi.org/10.1063/5.0161876
Google Scholar
Crossref
Search ADS
 
14.
H. G.
Lee
,
S.
Ham
, and
J.
Kim
, “
Isotropic finite difference discrettization of Laplacian operator
,”
Appl. Comput. Math.
22
(
2
),
259
–
274
(
2023
).
https://doi.org/10.30546/1683-6154.22.2.2023.259
Google Scholar
Crossref
Search ADS
 
15.
S.
Lee
,
S.
Yoon
, and
J.
Kim
, “
A linear convex splitting scheme for the Cahn–Hilliard equation with a high-order polynomial free energy
,”
Int. J. Numer. Methods Eng.
124
,
3586
–
3602
(
2023
).
https://doi.org/10.1002/nme.7288
Google Scholar
Crossref
Search ADS
 
16.
S.
Lee
,
S.
Yoon
,
C.
Lee
,
S.
Kim
,
H.
Kim
,
J.
Yang
,
S.
Kwak
,
Y.
Hwang
, and
J.
Kim
, “
Effective time step analysis for the Allen–Cahn equation with a high-order polynomial free energy
,”
Int. J. Numer. Methods Eng.
123
(
19
),
4726
–
4743
(
2022
).
https://doi.org/10.1002/nme.7053
Google Scholar
Crossref
Search ADS
 
17.
Y.
Li
 and
S.
Guo
, “
Triply periodic minimal surface using a modified Allen–Cahn equation
,”
Appl. Math. Comput.
295
,
84
–
94
(
2017
).
https://doi.org/10.1016/j.amc.2016.10.005
Google Scholar
Crossref
Search ADS
 
18.
Y.
Li
,
K.
Qin
,
Q.
Xia
, and
J.
Kim
, “
A second-order unconditionally stable method for the anisotropic dendritic crystal growth model with an orientation-field
,”
Appl. Numer. Math.
184
,
512
–
526
(
2023
).
https://doi.org/10.1016/j.apnum.2022.11.006
Google Scholar
Crossref
Search ADS
 
19.
Y.
Li
,
Q.
Xia
,
S.
Kang
,
S.
Kwak
, and
J.
Kim
, “
A practical algorithm for the design of multiple-sized porous scaffolds with triply periodic structures
,”
Math. Comput. Simul.
220
,
481
–
495
(
2024
).
https://doi.org/10.1016/j.matcom.2024.02.004
Google Scholar
Crossref
Search ADS
 
20.
C.
Liu
,
Z.
Qiao
, and
Q.
Zhang
, “
Multi-phase image segmentation by the Allen–Cahn Chan–Vese model
,”
Comput. Math. Appl.
141
,
207
–
220
(
2023
).
https://doi.org/10.1016/j.camwa.2022.12.020
Google Scholar
Crossref
Search ADS
 
21.
T.
Tang
 and
J.
Yang
, “
Implicit-explicit scheme for the Allen–Cahn equation preserves the maximum principle
,”
J. Comput. Math.
34
(
5
),
451
–
461
(
2016
).
https://doi.org/10.4202/jcm.1603-m2014-0017
Google Scholar
 
22.
B.
Xia
,
Y.
Li
, and
Z.
Li
, “
Second-order unconditionally stable direct methods for Allen–Cahn and conservative Allen–Cahn equations on surfaces
,”
Mathematics
8
,
1486
(
2020
).
https://doi.org/10.3390/math8091486
Google Scholar
Crossref
Search ADS
 
23.
B.
Xia
,
X.
Xi
,
R.
Yu
, and
P.
Zhang
, “
Unconditional energy-stable method for the Swift–Hohenberg equation over arbitrarily curved surfaces with second-order accuracy
,”
Appl. Numer. Math.
198
,
192
–
201
(
2024
).
https://doi.org/10.1016/j.apnum.2024.01.005
Google Scholar
Crossref
Search ADS
 
24.
Q.
Xia
,
Y.
Liu
,
J.
Kim
, and
Y.
Li
, “
Binary thermal fluids computation over arbitrary surfaces with second-order accuracy and unconditional energy stability based on phase-field model
,”
J. Comput. Appl. Math.
433
,
115319
(
2023
).
https://doi.org/10.1016/j.cam.2023.115319
Google Scholar
Crossref
Search ADS
 
25.
Q.
Xia
,
G.
Sun
,
Q.
Yu
,
J.
Kim
, and
Y.
Li
, “
Thermal-fluid topology optimization with unconditional energy stability and second-order accuracy via phase-field model
,”
Commun. Nonlinear Sci. Numer. Simul.
116
,
106782
(
2023
).
https://doi.org/10.1016/j.cnsns.2022.106782
Google Scholar
Crossref
Search ADS
 
26.
W.
Xie
,
J.
Feng
,
Q.
Xia
,
J.
Kim
, and
Y.
Li
, “
Design of the shell-infill structures using a phase field-based topology optimization method
,”
Comput. Methods Appl. Mech. Eng.
429
,
117138
(
2024
).
https://doi.org/10.1016/j.cma.2024.117138
Google Scholar
Crossref
Search ADS
 
27.
W.
Xie
,
Q.
Xia
,
Q.
Yu
, and
Y.
Li
, “
An effective phase field method for topology optimization without the curvature effects
,”
Comput. Math. Appl.
146
,
200
–
212
(
2023
).
https://doi.org/10.1016/j.camwa.2023.06.037
Google Scholar
Crossref
Search ADS
 
28.
X.
Xu
,
Y.
Hu
,
Y.
He
,
J.
Han
, and
J.
Zhu
, “
High-order analysis of lattice Boltzmann models for the conservative Allen–Cahn equation
,”
Comput. Math. Appl.
146
,
106
–
125
(
2023
).
https://doi.org/10.1016/j.camwa.2023.06.034
Google Scholar
Crossref
Search ADS
 
29.
H.
Zhang
,
X.
Qian
, and
S.
Song
, “
Third-order accurate, large time-stepping and maximum-principle-preserving schemes for the Allen–Cahn equation
,”
Numer. Algorithms
95
(
3
),
1213
–
1250
(
2024
).
https://doi.org/10.1007/s11075-023-01606-w
Google Scholar
Crossref
Search ADS
 
30.
H.
Zhang
,
X.
Qian
,
J.
Xia
, and
S.
Song
, “
Efficient inequality-preserving integrators for differential equations satisfying forward Euler conditions
,”
ESAIM: Math. Modell. Numer. Anal.
57
,
1619
–
1655
(
2023
).
https://doi.org/10.1051/m2an/2023029
Google Scholar
Crossref
Search ADS
 
31.
H.
Zhang
,
J.
Yan
,
X.
Qian
,
X.
Gu
, and
S.
Song
, “
On the preserving of the maximum principle and energy stability of high-order implicit-explicit Runge–Kutta schemes for the space-fractional Allen–Cahn equation
,”
Numer. Algorithms
88
,
1309
–
1336
(
2021
).
https://doi.org/10.1007/s11075-021-01077-x
Google Scholar
Crossref
Search ADS
 
32.
H.
Zhang
,
J.
Yan
,
X.
Qian
, and
S.
Song
, “
Numerical analysis and applications of explicit high order maximum principle preserving integrating factor Runge–Kutta schemes for Allen–Cahn equation
,”
Appl. Numer. Math.
161
,
372
–
390
(
2021
).
https://doi.org/10.1016/j.apnum.2020.11.022
Google Scholar
Crossref
Search ADS
 
© 2025 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).
This work is licensed under a Creative Commons Attribution 4.0 International License .
Attribution