Unconditionally stable monte carlo simulation for solving the multi-dimensional Allen–Cahn equation

1.   Introduction

In this paper, we propose an explicit hybrid numerical method for solving the multi-dimensional Allen–Cahn (AC) equation. The AC equation has emerged as a widely studied mathematical model, playing a pivotal role in the analysis of phase transitions ^[1]. This equation serves as a powerful tool for understanding the intricate dynamics that govern the boundary between ordered and disordered states in various physical systems ^[2]. With its ability to capture the evolution and behavior of interfaces, the AC equation has proven instrumental in fields, such as materials science, condensed matter physics and pattern formation ^[3,4]. Its innovative formulation has sparked numerous advancements in our comprehension of complex phenomena, shaping our understanding of phase transition dynamics and paving the way for further exploration and discovery in diverse scientific disciplines ^[5,6]. The multi-dimensional AC equation is given by

where ${\bf n}$ is outer unit normal vector and $\Omega \subset \mathbb{R}^d$ is a domain. Here, $d$ is the dimension of space, $\phi({\bf x}, t)$ is the concentration, $F(\phi) = 0.25\phi^2(\phi-1)^2$ and $\epsilon$ is an interface layer parameter. The AC equation is the $L^2$-gradient flow of the following total free energy functional:

In the AC Eq (1.1), we have the following diffusion equation:

which causes numerical difficulties when the space dimension is high, i.e., the curse of dimensionality which makes it computationally intractable because of lack of computer memory and excessive computational cost. Classical discretization schemes encounter a significant challenge known as the curse of dimensionality, which was initially introduced by Bellman ^[7] in the realm of optimal control problems for classical multidimensional partial differential equation (PDE) models ^[8]. The curse of dimensionality presents a pervasive challenge across diverse domains, such as mathematics, physics, data science and any context that deals with the arrangement, aggregation and examination of spaces with numerous dimensions ^[9,10,11]. This predicament becomes particularly conspicuous when endeavoring to tackle extensively utilized and extensively studied second-order semi-linear parabolic PDEs within high-dimensional spaces. There have been many studies to solve Eq (1.2) numerically. Shawn Koohy et al. ^[11] developed a method to mitigate the curse of dimensionality based on backward stochastic differential equations and deep neural networks to solve multi-dimensional AC equations. In ^[12], the parabolic sine-Gordon equation was solved using an operator splitting method that splits the linear and nonlinear terms, where the diffusion term is solved using the Fourier spectral method. Ayub et al. ^[13] solved the diffusion term in the AC equation using the Fourier spectral method. Lee ^[14] proposed the second-order operator splitting Fourier spectral method for the factional-in-space reaction-diffusion equation and solved the diffusion equation using the discrete cosine transform. Jeong et al. ^[15] proposed an explicit hybrid finite difference method using the explicit Euler method for the diffusion equation and a closed-form analytical solution for the nonlinear term by splitting the AC equation into a diffusion term and a nonlinear term based on the operator splitting method.

The stochastic heat equation is one of the representative equations among the stochastic partial differential equations. In ^[16], the authors discussed numerical positivity and almost surely exponential stability of the stochastic heat equation. Sun and Kumar ^[17] developed a numerical method for high-dimensional stationary Fokker–Planck equations using a tensor decomposition approach, focusing on the curse of dimension. Shrestha ^[18] solved numerically one-dimensional diffusion equation using the Monte Carlo simulation (MCS). The Fokker–Planck and Langevin equations with no drift term and constant diffusion coefficient are the diffusion equation. Medved et al. ^[19] compared and analyzed the Fokker–Planck and Langevin equations. Naeimi and Farshad ^[20] developed the finite volume Monte Carlo (FVMC) method for solving three-dimensional steady-state heat equations. In ^[21], Naeimi considered the application of MCS to a general form of heat equations and also investigated the floating random walk procedure to solve multi-dimensional problems with various boundary conditions. Nakagawa ^[22] verified a proposed method by presenting the results of numerical experiments applying the MCS to the stochastic heat equation with various boundary conditions. Venkiteswaran and Junk ^[23] developed a Quasi-Monte Carlo (QMC) method for diffusion equations in high dimensions and the numerical results showed faster convergence than the standard MCS. In ^[24], Novikov et al. explained Laplacian as a random for Brownian particle location and analyzed and presented the particle-based numerical method, MCS of Brownian diffusion effects, for transport equations for spheres (or circles).

Our proposed method is based on the operator splitting scheme. Various methods have been developed for solving the AC type equations based on the operator splitting method. In ^[25], Lee developed a numerical solution us ing the operator splitting scheme for the energy-dissipative and mass-conservative AC equation. The operator splitting method for the AC equation proposed in ^[26] is second-order accurate and stable. Cheng et al. ^[27] presented fast and stable operator splitting methods for phase-field models. Chertock et al. ^[28] proposed a fast explicit operator splitting method for the convection-diffusion equation. In addition, various studies on the AC equation have been performed by researchers. Poochinapan and Wongsaijai ^[29] developed a fourth-order compact structure-preserving difference method and verified that the numerical method has fourth-order and second-order accuracy in space and time, respectively. The primary purpose of this paper is to reduce the curse of dimensionality of the numerical algorithms for the AC equation. To mitigate the curse of dimensionality, our proposed method uses the MCS for the diffusion term and a closed-form analytic solution for the nonlinear term. Kai and Wei ^[30] developed high-order energy-stable methods for the multi-length-scale incommensurate phase-field crystal equation, which permits studying the phase behavior of aperiodic structures. The Cahn–Hilliard (CH) equation is one of the famous phase field models, along with the AC equation. Liupeng and Yunqing ^[31] developed a second-order scalar auxiliary variable approach in time and employed the linear finite element method in space to solve the CH type equation of the phase field crystal model.

The advantages of the proposed method are as follows:

● It can overcome the dimensional curse in solving the multi-dimensional AC equation.

● As the dimensionality increases, the convergence speed is faster compared to that of the finite difference method.

The contents of this paper are as follows. A detailed numerical solution algorithm is described in Section 2 and then the proposed method is verified through several numerical experiments in Section 3. Finally, a conclusion is presented in Section 4.

2.   Numerical algorithm

Using the operator splitting method, the proposed unconditionally stable MCS of the multi-dimensional AC equation consists of two steps. First, we solve the following diffusion equation using MCS.

Second, we solve the following ordinary differential equation (ODE) using a closed-form analytic solution ^[4].

Instead of solving the PDE (2.1), we use the following stochastic differential equation (SDE):

where ${\bf X}(t)$ is the position of a particle at time $t$, $d$ is the dimension of the space and $d{\bf W}(t) = \sqrt{dt} {\bf Z}$ is the Brownian motion. Here, ${\bf Z}$ is a random vector variable drawn from a standardized multivariate normal distribution. Then, Eq (2.1) is the Fokker–Planck equation of Eq (2.3) ^[18].

For simplicity of exposition, let us consider the solution algorithm in two-dimensional space, $\Omega = (L_x, R_x)\times(L_y, R_y)$. Let the spatial steps be $h_x$, $h_y$, the temporal step be $\Delta t$, the total number of particles be $M$ and $N_x = (R_x-L_x)/h_x$, $x_i = L_x+(i-0.5)h_x$ for $i = 1, \dots, N_x$, $N_y = (R_y-L_y)/h_y$, $y_j = L_y+(j-0.5)h_y$ for $j = 1, \dots, N_y$. We denote the numerical approximations of $\phi(x_i, y_j, n\Delta t)$ by $\phi_{ij}^{n}$, the position of $k$-th particle at time $n \Delta t$ by ${\bf X}_{k}^{n}$ and the value of the $k$-th particle at time $n\Delta t$ by $\psi_k^n$. The initial conditions for the particle positions ${\bf X}_k^0$, $k = 1, 2, \dots, M$ and the corresponding particle values $\psi_k^0$, $k = 1, 2, \dots, M$, are given as follows. If $\phi(x_i, y_j, 0) > 0.05$ for $i = 1, 2, \dots, N_x$, $j = 1, 2, \dots, N_y$, it is given as

where $q$ is some integer and $m$ is the number of particles given at per point. The particle position $X_{k}^n$ and particle value $\psi_k^n$ form pairs for $k = 1, 2, \dots, M$. For a better understanding, a schematic diagram showing the positions of some initial particles ${\bf X}_k^0$ for $k = 1, \dots, m$ in one-dimensional space is presented in Figure 1.

First, we solve the discretized diffusion equation. We update the particle positions ${\bf X}_k^n$ for $k = 1, \dots, M$ according to Eq (2.3) as follows:

Second, solve for the discretized nonlinear Eq (2.2). We define the local sum of particle values as follows:

where $\Omega_{ij} = (L_x+h_x(i-1), L_x+h_xi)\times (L_y+h_y(j-1), L_y+h_yj)$ is a cell for some $i$ and $j$. Figure 2(d), (e) and (g) schematically show this procedure. If $\phi_{ij}^* > tol$, $i = 1, 2, \dots, N_x$, $j = 1, 2, \dots, N_y$, then we analytically solve the following ODE (2.2):

with an initial condition $\Phi(0) = \phi_{ij}^*$, where $tol$ is a given tolerance and is used to avoid division by an extremely small number. The analytic solution of Eq (2.6) with the initial condition $\phi_{ij}^*$ after time $\Delta t$ is given as

Then, we update the particle value $\psi_k^n$ for the $k$-th particle position ${\bf X}_k^n$ for $k = 1, 2, \dots, M$ as follows:

This means that we have analytically solved Eq (2.6) with the initial condition $\Phi(0) = \phi_{ij}^*$ for each cell $\Omega_{ij}$ as follows:

Figure 2(f), (h) show the results of the final procedure when the value of $\phi_{ij}^*$ is less than $0.5$ or greater than $0.5$, respectively.

We note that once $\psi_k^{n+1} > 0.5$, then $\psi_k^{\hat n} > 0.5$ is satisfied for all time, ${\hat n} > n+1$. To avoid this situation, if $\psi_k^{n+1} > 0.5$, then we add one additional particle ${\bf X}_{M+1}^{n+1}$ at the same position ${\bf X}_k^{n+1}$, i.e., ${\bf X}_{M+1}^{n+1} = {\bf X}_k^{n+1}$. In addition, we divide the value of particle $\psi_k^{n+1}$ by $2$, i.e., $\psi_k^{n+1} = \psi_k^{n+1}/2$ and set $\psi_{M+1}^{n+1} = \psi_k^{n+1}$.

3.   Numerical experiments

For an appropriate interfacial parameter value, we consider an equilibrium solution of the AC equation, which is given by $\phi(x) = 0.5(1+\tanh(x/(2\sqrt{2}\epsilon)))$. Based on this equilibrium solution, we define the interfacial layer parameter $\epsilon_l$ as follows:

which implies that we have approximately $lh$ width across the interfacial transition layer ^[32]. The schematic diagram for the interfacial layer is presented in Figure 3.

3.1.   Stability tests

In this section, we demonstrate the stability of the proposed algorithm. For simplicity of proof, we consider the proposed algorithm in two-dimensional space as in Section 2. Similarly, we can prove this in another multi-dimensional space. We verify the unconditional stability and energy decreasing for the proposed algorithm through numerical experiments. We have used the zero Neumann boundary condition as

The discrete total energy ${\mathcal{E}}_h^d(\phi^n)$, discrete total mass ${\mathcal{M}}_h^d(\phi^n)$ and discrete maximum value $Max^d(\phi^n)$ are defined as:

The first test is a stability test of the proposed method for the spatial steps. In this section, we use the following initial condition on a computational domain $\Omega = (-3, 3)\times(-3, 3)$

Here, the parameters used are $\Delta t = 0.0025$, $T = 0.5$, $\epsilon = 0.2$ and different spatial steps $h = 0.2, \; 0.1, \; 0.05$ with corresponding $m = 1000, \; 2000, \; 4000$, respectively. Numerical experiments were performed to assess the stability for different spatial steps $h = 0.2, \; 0.1, \; 0.05$ and the results of the stability test are depicted in Figure 4.

As a second test, we consider the stability of the temporal step. We used parameters $h = 0.05$, $T = 0.5$, $m = 1400$, $\epsilon = \epsilon_{10}$ and different temporal steps $\Delta t = 2h^2, \; h^2, \; 0.5h^2$. Figure 5 (a)–(c) show the results of stability tests as discrete total energy, discrete total mass and discrete maximum value, respectively.

We numerically performed spatial and temporal stability tests and observed that the numerical solution does not explode over various spatial and temporal steps.

3.2.   Convergence test

The proposed method is affected by the initial total number of particles $M$. However, the $M$ required to obtain a numerical solution up to the total time $T$ is theoretically impossible or difficult to obtain. Therefore, we performed a numerical test for the initial number of particles per point $m$ required to obtain the enough $M$ in two-and three-dimensional spaces. Note that this numerical test is a test to show the approximate $m$ required for the proposed method. The numerical test repeats the numerical simulation by adding $5$ to $m$ until the absolute error between the analytical solution and the numerical solution of the proposed method is less than tolerance $0.001$. The parameters used are $r = 0.8$, spatial step $h = 0.05$, total time $T = 0.01$, $\epsilon = \epsilon_8$ and $x_i = -2.025+0.05j$ for $i = 1, \dots, 80$, $y_j = -2.025+0.05j$ for $j = 1, \dots, 80$. Here, $m$ particles are given to the $(x_i, y_j)$ when $\phi(x_i, y_j, 0)$ is greater than $0.05$ at the $(x_i, y_j)$, therefore the initial total number of particles $M$ becomes 1264 $m$ and 33,552 $m$ in two-and three-dimensional space, respectively. Figure 6(a), (b) show the absolute errors of numerical and analytic solutions in two-and three-dimensional space for the initial number of particles per point $m$, respectively. As a result of the numerical test, the m required for the absolute error to be smaller than the absolute error is $200$ in two-dimensional space and $105$ in three-dimensional spaces. We observe that the absolute error decreases sharply as $m$ increases and the required $m$ decreases as the dimension of space increases. However, because the initial total number of particles $M$ used is $252,800$ in the two-dimensional space and $3,522,960$ in the three-dimensional space, it increases as the dimension of space increases.

3.3.   Motion by mean curvature and ratio of CPU times

As $\epsilon$ goes to zero, the zero level set of the solution of the AC equation converges to the mean curvature flow ^[33,34]. Thus, we consider the mean curvature flow of $(d-1)$-dimensional spheres in $d$-dimensional space. Let us define the gamma function as

thus $\Gamma(k) = (k-1)!$, $\Gamma(k+1/2) = (k-1/2)(k-3/2)\cdots(1/2)\sqrt{\pi}$, for a positive integer $k$. Then, the area of sphere with radius $R = 1$ at the origin is $A(R) = 2\pi^{\frac{d}{2}} / \Gamma \left(\frac{d}{2}\right).$ We define a $d$-ball as a region comprising of all points within a radius $R$ of a given point in $d$-dimensional space. The volume of $d$-ball with radius $R = 1$ at the origin is

Therefore, we obtain volume of $d$-ball with any radius $R > 0$ at the point ${\bf x}\in{\mathbb{R}^d}$ as $V(R) = 2\pi^{\frac{d}{2}} R^d / [d\, \Gamma(\frac{d}{2})].$ More details can be found in ^[35]. Then, we can express the radius $R$ in terms of $V$:

In $d$-dimensional space, the radius of $d$-ball can be calculated using the volume of $d$-ball and the calculated radius becomes the radius of the $(d-1)$-dimensional sphere. Then, mean curvature of sphere of radius $R$ is $H = (d-1)/R$. The mean curvature flow equation reduces to the following ordinary differential equation, for a sphere of initial radius $R_0$:

The solution of Eq (3.2) is $R(t) = \sqrt{R^2_0-2(d-1)t}$. The volume of the numerical solution is obtained as follows: $V(n\Delta t) = \sum_{k = 1}^M \psi_k^n h^d.$ Therefore, by substituting the obtained volume $V(n\Delta t)$ into Eq (3.1), the radius $R(n\Delta t)$ is obtain. Because the proposed method is used in high-dimensional, instead of calculating the radius using the distance from the center to the zero-contour, we use the volume of the solution to calculate the radius. In numerical experiments, the initial radius $R_0$ is calculated numerically using the volume of the initial condition. Consequently, we assumed that for the analytic solution of the zero level set of the AC equation for the motion by mean curvature, $\epsilon$ becomes zero.

3.3.1.   Two-dimensional space

The initial condition is given as

on the computational domain $\Omega = (-2, 2)\times(-2, 2)$, where $r$ is radius. The analytic solution of the zero level set of the AC equation, as $\epsilon$ converges to zero, is given by:

which represents the analytic solution for the motion by mean curvature of a $1$-dimensional sphere in $2$-dimensional space. The parameters used are $r = 0.8$, spatial step $h = 0.05$, total time $T = 0.1$, $\epsilon = \epsilon_8$, $m = 600$ and $x_i = -2.025+0.05j$ for $i = 1, \dots, 80$, $y_j = -2.025+0.05j$ for $j = 1, \dots, 80$. Because the initial radius is $0.8$, the analytic solution to the motion by mean curvature of the AC equation is $R(t) = \sqrt{0.64-2t}$. To verify the temporal convergence of the proposed method and compare the explicit Euler method, we used different temporal steps $\Delta t = 10h^2, \; 5h^2, \; 2.5h^2$ for the proposed method and $\Delta t = 0.1h^2$ for the explicit Euler method. We used the initial number of particles per point $m = 600$, thus the initial particles are given $600$ particles at position $(x_i, y_j)$ if $\phi(x_i, y_j, 0)$ is greater than $0.05$ at the cell center $(x_i, y_j)$. Therefore, the initial number of particles $M = 758,400$ for numerical experiments in two-dimensional space. In addition, for explicit Euler method, a value was given only when $\phi(x_i, y_j, 0)$ at the cell center $(x_i, y_j)$ was greater than $0.05$ to make the initial radius $R_0$ the same. We consider the effects of temporal steps on the dynamics of the AC equation. Figure 7 shows temporal evolutions of the radius $R(t)$ using the proposed method with different temporal steps $\Delta t = 10h^2, \; 5h^2, \; 2.5h^2$ and the explicit Euler method with temporal step $\Delta t = 0.1h^2$ in two-dimensional space. As the temporal step $\Delta t$ decreases, we can observe that the numerical solutions converge to the analytic solution $R(t)$. In addition, it can be seen that the proposed method is unconditional stable with respect to time when a enough initial total number of particles $M$ is given.

To verify the computational speed of the numerical algorithm, we defined CPU time as the elapsed time of the main algorithm and defined the CPU time$_{Ratio}$ as the value obtained by dividing the CPU time of the numerical method by the CPU time of the explicit Euler method. Table 1 lists the CPU time and CPU time$_{Ratio}$ to solve the AC equation up to total time $T = 0.1$ for the proposed method with different temporal steps and the explicit Euler method in two-dimensional space. In low dimensional space, such as two-dimension space, the explicit Euler method computes faster than the proposed method. However, we perform numerical experiments in three- and four-dimensional space to show that our proposed method is more efficient in solving the AC equation as the dimension of space increases.

Unless otherwise specified, parameters used in two-dimensional space expand equally in three- and four-dimensional spaces.

3.3.2.   Three-dimensional space

Numerical experiments similar to those in two-dimensional space are performed in three-dimensional space as well. We consider the initial conditions as

on computational domain $\Omega = (-2, 2)\times(-2, 2)\times(-2, 2)$. The analytic solution of the zero level set of the AC equation, as $\epsilon$ converges to zero, is given by:

which represents the analytic solution for the motion by mean curvature of a $2$-dimensional sphere in $3$-dimensional space. The parameters used in the three-dimensional space are extensions of the parameters used in two-dimensional except for the initial number of particles per point $m = 300$ and therefore initial number of particles $M = 10,065,600$. Because the initial radius is $0.8$, the analytic solution to the motion by mean curvature of the AC equation is $R(t) = \sqrt{0.64-4t}$.

Figure 8 shows the temporal evolution of radius of the proposed method with difference temporal steps $\Delta t = 10h^2, \; 5h^2, \; 2.5h^2$ and explicit Euler method with temporal step $\Delta t = 0.1h^2$ in three-dimensional space. We can observe that, given an enough initial total number of particles $M$ in three-dimensional space, the proposed method generates good results despite of the large temporal steps $\Delta t$.

Table 2 lists the CPU time and CPU time$_{Ratio}$ to solve the AC equation up to total time $T = 0.1$ for the proposed method with different temporal steps and the explicit Euler method in three-dimensional space. Compared with Table 1, we can observe that the overall CPU time$_{Ratio}$ is reduced, which means that the efficiency of the proposed method in three-dimensional space is better than that in two-dimensional space.

3.3.3.   Four-dimensional space

We also perform the same numerical experiments on the four-dimensional space. The initial condition is given as

in computational domain $\Omega = (-2, 2) \times(-2, 2)\times(-2, 2)\times(-2, 2)$. The analytic solution of the zero level set of the AC equation, as $\epsilon$ converges to zero, is given by:

which represents the analytic solution for the motion by mean curvature of a $3$-dimensional sphere in $4$-dimensional space. The parameters used in four-dimensional space are the same to two-and three-dimensional spaces except for the initial number of particles per point $m = 150$ and therefore initial number of particles $M = 118,560,000$. Because the initial radius is $0.8$, the analytic solution to the motion by mean curvature of the AC equation is $R(t) = \sqrt{0.64-6t}$.

In Figure 9, we can observe that the numerical solution of the proposed method converges to the analytical solution when the temporal step $\Delta t$ is decreased.

Table 3 lists the CPU time to solve the AC equation up to total time $T = 0.1$ using the proposed method with different time steps and for the explicit Euler method in four-dimensional space and the CPU time$_{Ratio}$ accordingly. Contrary to the numerical experiment on the computational speed in two-dimensional space, we can observe that the CPU time of the proposed method is faster than the CPU time of the explicit Euler method.

Next, we perform numerical experiments to demonstrate that the proposed method can solve the AC equation even for a long time, provided that the total number of particles is sufficient. The initial condition is given as

in computational domain $\Omega = (-3, 3) \times(-3, 3)\times(-3, 3)\times(-3, 3)$. The parameters used are $r = 2$, $T = 0.6$, $h = 0.1$, $\Delta t = 5h^2$, $\epsilon = \epsilon_{10}$ and $m = 70$. Consequently, the total initial number of particles was $M = 134,482,880$ according to the initial number of particles per point ($m$). Figure 10 shows the temporal evolution of the radius using the proposed method and the explicit Euler method for a long time.

We observed that the proposed method solves the AC equation well when an adequate initial number of particles is provided.

3.4.   Other nonlinear term

We consider the following Fisher–Kolmogorov–Petrovsky–Piscounov (Fisher–KPP) equation in two-dimensional space with $\phi(1-\phi)$ as an other nonlinear term ^[36,37].

We modify Eq (2.6) of the numerical algorithm to solve the nonlinear term of the Fisher–KPP equation using the proposed method.

with an initial condition $\Phi(0) = \phi_{ij}^*$. Equation (3.4) is a logistic growth model and an analytic solution can be obtained ^[38]. The analytic solution of Eq (3.4) with the initial condition $\phi_{ij}^*$ after time $\Delta t$ is given as

To solve Eq (3.3) in the domain $(-20, 20) \times (-20, 20)$, we set the computational domain to $(-25, 25)\times(-25, 25)$ because the boundary must be considered for the initial condition.

The exact solution is given by

For detailed information on the exact solution, please refer to ^[39]. The parameters used are $h = 2$, $T = 4$, $\Delta t = 0.01h^2$ and $m = 10,000$. We conduct numerical experiments by adjusting the standard value of $0.05$, which provides particles at the beginning of the proposed algorithm, to $0.001$ while considering smooth initial conditions. Figure 11 shows the temporal evolution of the numerical solution of the Fisher–KPP equation.

We observe that the proposed method can be applied to an other nonlinear equation with the nonlinear term which is different from $-F'(\phi)/\epsilon^2$.

4.   Conclusions

In this paper, we presented an efficient and novel unconditionally stable MCS for solving the multi-dimensional AC equation that can model the motion by mean curvature flow of a hypersurface. We used an operator splitting method in which the diffusion and nonlinear terms are separately solved. The diffusion term is calculated using MCS for the stochastic differential equation and the nonlinear term is locally computed for each particle in the virtual grid. Therefore, the proposed method is unconditionally stable with respect to time, given an enough number of particles to solve the AC equation up to the total time. We performed numerical experiments in two-, three- and four-dimensional spaces to verify the performance of the proposed algorithm. The proposed algorithm is more efficient in high-dimensional space because the algorithm can be easily scaled despite the increase in the space dimension. Through numerical experiments, we observed that the proposed algorithm solves the AC equation more efficiently as the dimension of space increases.

Use of AI tools declaration

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

Acknowledgments

The corresponding author (J. S. Kim) expresses thanks for the support from the BK21 FOUR program. The authors express their gratitude to the reviewers for their valuable and insightful feedback on the revised version of this article.

Conflict of interest

The authors declare there are no conflicts of interest.