Computing the area-minimizing surface by the Allen-Cahn equation with the fixed boundary

1.   Introduction

The Allen–Cahn (AC) equation has its origin in material science ^[1], but it has been studied in many fields such as partial differential equations ^[3], geometry ^[8], scientific computations ^{[4,10,13,19,26,27]}, image processing ^[17], computational biology ^[2,23] and references therein. Introducing the order parameter $\phi$ and positive and small parameter $\epsilon$ in the AC equation, we can observe the constant functions $\pm\phi$ as the solution. That is, we have the domain $\Omega = \Omega_{+}\cup \mathcal{I} \cup \Omega_{-}$, where $\Omega_{+}: = \{\mathbf{x}\in \mathbb R^3|\phi(\mathbf{x}) = 1\}$, $\Omega_{-}: = \{\mathbf{x}\in \mathbb R^3|\phi(\mathbf{x}) = -1\}$, and $\mathcal{I}: = \{\mathbf{x}\in \mathbb R^3|-1 < \phi(\mathbf{x}) < 1\}$. It is of interest to analyze the configuration when two phases coexist. The two phases are distinct from each other by taking one of the functions $\phi = \pm1$. The constant function takes values close to $1$ in the subdomain $\Omega_{+}$, and the other value takes $-1$ in the subdomain $\Omega_{-}$ as illustrated in Figure 1(a). One of the reasons that the AC equation is famous for many applications that it has the intrinsic behavior of phase separation $\pm 1$ in the domain $\Omega$.

Phase separation is the process of a single homogeneous mixture of two or more components spontaneously dividing into two or more distinct phases. The components of a mixture can be separated into different phases based on their different properties, e.g., their solubility or density. The process can occur in a variety of systems, including fluid dynamics ^[29,30].

Mathematically, the phase separation induced by the AC equation can be interpreted as the gradient flow of the Ginzburg-Landau energy functional. For the gradient flow which occurs with phase separations, the behavior of the interface between two phases follows the motion by curvature as the interfacial width $\epsilon$ goes to $0$. As the singular limit, the curvature driven motion provides a link between the solution from the AC equation and minimal hypersurface. Thus, the AC equation has been scrutinized analytically and numerically for geometric analysis ^[12,15].

The Allen-Cahn (AC) equation is a semilinear PDE which is closely related to the theory of minimal surface ^[8]. The AC equation is as follows:

Here, $\phi$ is an order parameter, $F(\phi)$ is the double-well free energy as shown in Figure 1(b), where we have two minima $\phi = -1$ and $\phi = 1$, and $F'(\phi)$ is the derivative of $F(\phi)$ with respect to $\phi$. The double-well potential $F$ and its derivative with respect to $\phi$ are given by

In terms of the frame constraints to surface, the free angle is implemented at the frame boundaries that surface is attached to the fixed frame all the time, but the angle between surface and frame is free. We implicitly prescribe the frame boundary conditions. but angles in the boundaries are free by maintaining the zero interface of $\phi$ in the frame boundaries.

Constructing a minimal surface via a numerical method has been studied for decades ^[6,9,14]. Most numerical schemes use the level-set. Their works manage to constrain a surface to the frame boundaries by reattaching the surface to the frame boundary iteratively or enforce the surface sticking to the frame boundaries by imposing zero values repeatedly. Yet, in other parts of numerical methods, one can construct minimal surfaces as the weak limit of level sets of a semilinear elliptic PDE. Thus, one of the well-known equations is the AC equation. To secure the frame boundary conditions, we combine the AC's Laplacian operator with the constraint operator. In this way, we can avoid the instability in the frame boundaries, instead of repeatedly imposing the zero values in the frame boundaries explicitly. On the other hand, there are works related the AC equation and minimal surface ^[16,18,31]. We notice that their numerical works only concern minimal surfaces with free boundary conditions, e.g., triply periodic minimal surfaces.

It is worth mentioning that minimal surfaces find extensive utility across numerous industrial applications. In the research work ^[20], the composite scaffold can be formulated by integrating two distinct triply periodic minimal surfaces. Along with advances in additive manufacturing, porous structures are easily manufactured. In that sense, the modified phase-field equation is a popular choice for multiscale and minimizing area topology optimization of porous structures ^[32].

In the differential geometric aspect, minimal surfaces as soap films are beautiful geometric objects that can be observed everywhere in nature. Researchers have studied them for the past years with the help of mathematical theories, experiments and even recently computer simulations. Soap film have many interesting properties such as geometric shapes and macroscopic and molecular behaviors. It is minimal because the surface tension reduces its area to a minimum. For that reason, mathematical background for soap films with the frame is mostly related to minimal surface theory. Minimal surfaces have their root in the calculus of variations developed by Euler and Lagrange and in later investigations that have been performed over a few centuries.

Euler discovered a surface when a catenary is rotated around the $x$-axis, and one gets a surface which minimizes surface area. This surface is called the catenoid. Later, Meusnier formulated it as a solution of Lagrange's equation. The catenoid has genus zero, two ends and total curvature $4\pi$, i.e., the only properly embedded minimal surface with two ends and finite curvature is the catenoid. Embedded minimal surfaces with ends have been studied in geometry ^[5,21,22]. A catenoid is created via the equilibrium shape of a soap film, which is stretched between two parallel circular rings. On the one hand, Scherk in 1834 discovered a doubly periodic, complete, and properly embedded minimal surface in Euclidean space. This surface is non-trivial and also spanned into a frame. In fact, Scherk's surface is either parameterized by punctured spheres and then has one translational period or one screw motion period or it is also parameterized by rectangular tori, implying being doubly periodic. Catenoid and Scherk's surface are illustrated in Figure 2.

The main contribution of this work is to stably secure the frame boundaries so that the area-minimizing surface with the fixed frame boundary is easily and more efficiently constructed. We believe that the numerical treatment for the frame boundaries provides a foundation to proceed to much deeper understanding of minimal surfaces and the behavior of the AC equation.

The rest of this manuscript is organized as follows. In Section 2, we present a brief review of differential geometry about the Weierstrass-Enneper representation. By using the representation, we can compare the analytic surface with the numerical surface, and also verify that our numerical surface is minimal. In section 3, we provide the numerical treatment for the frame boundary conditions, and then perform numerical simulations to validate the proposed scheme. In the conclusion, section 5, we summarize our work with further discussions.

2.   Theoretical approach for minimal surface

In this paper, we start with the differential geometry knowledge that can be used for numerical verification of the proposed scheme. We give the basics of the Weierstrass-Enneper integral representation of minimal surfaces, demonstrating that Scherk's first surface and catenoid can be described in terms of holomorphic functions of a complex variable. Later, we compare the solutions from the Weierstrass representation with numerical solutions in order to verify our proposed numerical scheme.

We consider that $\mathbf X:U\to \mathbb R^3$ is a parameterized surface of the form:

With the parameterization $\mathbf X(u, v)$, it is said to be isothermal if

A minimal surface is a surface $\mathcal{M}$ with mean curvature $H\equiv 0$ at all points $p\in \mathcal{M}$. When we suppose $E = \mathbf X_u\cdot \mathbf X_u$, $F = \mathbf X_u\cdot \mathbf X_v$, $G = \mathbf X_v\cdot \mathbf X_v$, as the coefficients of the first fundamental form, and $l = \mathbf X_{uu}\cdot \textbf{n}$, $m = \mathbf X_{uv}\cdot \textbf{n}$, and $n = \mathbf X_{vv}\cdot \textbf{n}$ as the coefficients of the second fundamental form. Then, the formula of mean curvature is $H = (lG-2mF+nE)/(2EG-2F^{2})$. If the parameterization $\mathbf X$ is isothermal and minimal, then we have

The Weierstrass-Enneper representation is well known for the useful way that it creates minimal surfaces. In order to know how it is possible, let $\mathcal{M}$ be a minimal surface defined by an isothermal parametrization $\mathbf X(u, v)$ and $z = u+iv$, corresponding complex coordinate, with $u = (z+\bar z)/{2}$ and $v = (-i \left(z-\bar z\right))/{2}$. So, we can write

Then, we can also know ${\mathbf{x}}(z, \bar z) = (x^{1}(z, \bar z), x^{2}(z, \bar z), x^{3}(z, \bar z))$, and $x^{i}(z, \bar z)$ is a complex valued function which takes real values. So, $x^{i}$ is ${\partial x^{i}}/{\partial z} = (x^{i}_{u}-ix^{i}_{v})/2$ by the definition of ${\partial}/{\partial z}$, and let us define

Since $(x^{i}_{z})^{2} = (1/4)((x^{i}_{u})^{2}-(x^{i}_{v})^{2}-2ix^{i}_{u}x^{i}_{v})$, and ${\mathbf{x}}(u, v)$ is isothermal,

As we know, if $(\phi)^{2} = 0$, the parametrization is isothermal. If we suppose $\mathcal{M}$ is a surface with parametrization ${\mathbf{x}}$, $\phi = \partial {\mathbf{x}}/\partial z$ and $(\phi)^{2} = 0$, $\mathcal{M}$ is minimal, and $\phi^{2}$ is holomorphic. In addition, $(\phi)^{2} = 0$ if and only if each $\phi^{i}$ is holomorphic. This implies that

Let $(\phi)^2 = 0$ and with the definition above $\phi = \partial \mathbf X/\partial z$, we have the following Weierstrass-Enneper representation.

Definition 1. The Weierstrass-Enneper representation I

Suppose $f$ is holomorphic on a domain $D$, $g$ is meromorphic on $D$, and $fg^{2}$ is holomorphic on $D$. Then, we can get a minimal surface defined by the parametrization ${\mathbf{x}}(z, \bar{z}) = (x^{1}(z, \bar{z}), x^{2}(z, \bar{z}), x^{3}(z, \bar{z}))$, where

Let us choose one function instead of two, and then we get a holomorphic $g$ which has an inverse function $g^{-1}$(that is holomorphic). So, we define $\tau = g$, which is $\, \mathrm{d}\tau = g' \, \mathrm{d} z$, and $F(\tau) = f/g'$, which is $F(\tau) \, \mathrm{d}\tau = f \, \mathrm{d} z$. Then, we can replace $g$ with $\tau$ and $f \, \mathrm{d} z$ with $F(\tau) \, \mathrm{d}\tau$.

Definition 2. For any analytic function $F(\tau)$, the Weierstrass-Enneper representation II is given by

where we define ${\mathbf{x}}(z, \bar{z}) = (x^{1}(z, \bar{z}), x^{2}(z, \bar{z}), x^{3}(z, \bar{z}))$.

According to the Weierstrass-Enneper representation II, we get that

Let us define $\tau = :u+i v$, and then we use $1+i\tau = 1-v+iu$ and $1-i\tau = 1+v-iu$. With $\tau = :u+i v$, the Weierstrass-Enneper representation II of Scherk's first surface is given by $F(\tau) = {2}/(1-\tau^4)$. Then, we can find the surface $({\mathbf{x}}_1(u, v), {\mathbf{x}}_2(u, v), {\mathbf{x}}_3(u, v))$:

Note that ${\mathbf{x}}_3$ can be computed in a similar fashion as in the part ${\mathbf{x}}_2$, and ${\mathbf{x}}_1$ can be rewritten by ${\mathbf{x}}_1 = \ln(\cos {\mathbf{x}}_3/\cos {\mathbf{x}}_2)$. Here, a piece of the surface is constructed within the square domain $\Omega: = ({\mathbf{x}}_2, \; {\mathbf{x}}_3)\in[-\pi/2, \pi/2]\times[-\pi/2, \pi/2]$.

Next, we see the next example. The Weierstrass-Enneper representation II of the catenoid is given by $F(\tau) = {1}/(2\tau^2)$. To see the surface, let $F(\tau) = {1}/(2\tau^2)$ be a holomorphic function. Then, with the substitution $\tau$ = $e^z$, we obtain the catenoid ${\mathbf{x}}(u, v)$. Together with holomorphic function $F(\tau)$, we get the catenoid by the parameterization ${\mathbf{x}}(u, v) = ({\mathbf{x}}_1(u, v), {\mathbf{x}}_2(u, v), {\mathbf{x}}_3(u, v))$:

For more details, we refer the reader to ^[7,24].

3.   Numerical implementation of the Allen-Cahn equation with the frame boundary

The finite difference method is one way of solving a differential equation using an approximate value of a derivative. In other words, this method is an approximate solution for $\phi({\mathbf{x}}, t)$ at a finite set of ${\mathbf{x}}$ and $t$. Note that we discretize the fixed boundary conditions in addition to the operator of the AC equation. Let us consider the $3$ dimensional domain $\Omega = (0, 1) \times (0, 1) \times (0, 1)$, and then discretize it as follows: $\Omega_h = \{ (x_i, \; y_j, \; z_k) | x_i = (i-0.5)h, \; y_j = (j-0.5)h, \; z_k = (k-0.5)h, \; i = 1, \ldots, N_x, \quad j = 1, \ldots, N_y \mbox{ and } k = 1, \ldots, N_z \}$, where $N = N_x = N_y = N_z$ is a large integer. Therefore, the computational domain $\Omega_h$ is defined by

Here, $h = 1/N = 1/N_x = 1/N_y = 1/N_z$ is the spatial step size. Let $\phi_{ijk}^n = \phi(x_i, y_j, z_k, t^n)$, where $t^n$ is the discrete time. Let $\Delta t = t^{n+1}-t^n$ be a temporal step. For simplicity of exposition, we denote $\phi({\mathbf{x}}_{i, j, k}, n\Delta t)$ by $\phi_{i, j, k}^n$, where $\Delta t = T/n_T$ is a temporal step size, $T$ is the final time, and $n_T$ is the total number of temporal steps. With these notations, let us consider the fully explicit Euler scheme for the AC equation with the time-dependent interfacial parameter:

where the three-dimensional discretization for the Laplace operator is given by

Here, the operator $\Delta_n$ indicates the Laplacian with the Neumann boundary condition. For the sake of notation convenience, we write the discrete operator as the the Neumann Laplacian. We let the interfacial width parameter $0 < \epsilon \ll 1$, and $\Omega\subset \mathbb R^d$ for $d = 3$, where $d$ is the number of dimensions. The other cases $d = 1, 2$ can be considered similarly in general. For constraints, we make the zero-level set of $\phi$ lie on $\Gamma\subset\Omega$. In fact, the AC equation can be interpreted as the gradient flow of the Ginzburg-Landau energy functional with constraints

Note that $\gamma > 0$, and $\mathcal{C}_g$ is dependent on the corresponding geometric constraints, which must satisfy the following condition:

This operator $C_g$ is for the frame boundary conditions. Its basis is to attach the surface to the frame by the boundary condition of the surface. Here, the operator $C_g$ is implicitly treated to improve the stability by combining the frame boundary condition with the discrete Laplacian operator. We elaborate on the discretized form of $\mathcal{C}_g$ in the subsequent paragraphs.

To solve Eq (3.1) with the fixed frame boundaries, the operator splitting method is applied:

If we vectorize $\phi_{i, j, k}^n$ for numerical purposes, the discrete solution at the $n$-th temporal step is simply represented by

First, we discretize Eq (3.2) as below:

We set up a symmetric matrix $C$ such that

where we define ${\mathbf{x}}_{i, j, k, p, q, r}: = (1-\omega_{p, q, r}){ {\mathbf{x}}_{i, j, k}+\omega_{p, q, r} {\mathbf{x}}_{p, q, r}}$. The positive coefficients $\alpha_{p, q, r} > 0$, $\beta_{p, q, r} > 0$ and $1 > \omega_{p, q, r} > 0$ are properly chosen to satisfy $\phi({\mathbf{x}}_{i, j, k, p, q, r}) = 0$ and symmetry of matrix $C$ elements.

Second, we solve equation (3.3) with $\phi^{*}$ analytically using the method of separation of variables introduced in ^[19,28], but we apply this for the $3$ dimensional case:

For the frame boundaries, we make use of the fixed boundary condition implicitly. Before further proceeding, we observe the Laplacian operator with the frame boundary condition as follows. We combine the geometric constraints with the AC's Laplacian operator for the fixed frame boundaries. Note that the discrete operator $\mathcal{L}_n$ approximates the Neumann Laplacian $\Delta_n$, and we define $C$ by the constraint operator for (3.5). To solve the system of Eqs (3.4) and (3.5), we combine the Neumann Laplacian $\mathcal{L}_n$ and boundary constraints $C$, i.e., $-\mathcal{L}_c: = -\mathcal{L}_n+C$. After vectorizing our computations, we obtain for $\phi^*$

Lemma 3. The Laplacian with constraints, $-\mathcal{L}_c$ is symmetric and positive semidefinite. Moreover, the linear operator $(I_n-\Delta t \mathcal{L}_c)$ is symmetric and positive definite.

Proof. As stated above, we describe the symmetric matrix $C$ such that

where ${\mathbf{x}}_{i, j, k, p, q, r}: = (1-\omega_{p, q, r}){ {\mathbf{x}}_{i, j, k}+\omega_{p, q, r} {\mathbf{x}}_{p, q, r}}$. The positive coefficients $\alpha_{p, q, r} > 0$, $\beta_{p, q, r} > 0$ and $1 > \omega_{p, q, r} > 0$ are properly chosen to satisfy $\phi({\mathbf{x}}_{i, j, k, p, q, r}) = 0$ and symmetry of matrix $C$ elements. Once the diagonal coefficient ${\bar{{ { \alpha_{i, j, k}} }}}$ is determined for the diagonal index $(i, j, k)$ of $C$, the off-diagonal coefficient ${\bar{{ { \beta_{p, q, r}} }}}$ follows $\left({\bar{{ { \beta_{i, j, k}} }}}^2/ {\bar{{ { \alpha_{p, q, r}} }}}\right) \phi_{p, q, r} = {\bar{{ { \beta_{i, j, k}} }}}\phi^*_{i, j, k}$ for the other diagonal index $(p, q, r)$ and off-diagonal index $(i, j, k)$ of $C$. This implies that the matrix $C$ is symmetric and singular, but the geometric constraints solely influence the specific elements of the matrix $C$. Since $-\mathcal{L}_c: = -\mathcal{L}_n+C$, we decompose it into the Neumann Laplacian $\mathcal{L}_n$ and frame boundary constraints $C$.

Let us first look at the eigenvalues of the Neumann Laplacian $\mathcal{L}_n$. With notations $\theta = \pi k/N_x$, $\psi = \pi l/N_y$, and $\rho = \pi l/N_z$, the eigenvalues of $-\mathcal{L}_n$ are of the form

Then, the combined linear operator $-\mathcal{L}_c$ is strictly diagonally dominant. Since a symmetric and diagonally dominant matrix $\left|(\mathcal{L}_c)_{i, i}\right| > \sum_{i\ne j} \left|(\mathcal{L}_c)_{i, j}\right|$ with real non-negative diagonal entries is positive definite, the Laplacian with constraints operator $-\mathcal{L}_c = -\mathcal{L}_n+C$ is also symmetric and positive semidefinite. Furthermore, $(I_n-\Delta t \mathcal{L}_c)$ is symmetric and positive definite for $\Delta t > 0$. □

Theorem 4. For $0 < \Delta t \leq h^2/(2d+2)$, the numerical scheme from (3.4), (3.5) and (3.6) preserve the boundedness of numerical solutions by $1$. In other words, $\max|\phi_{i, j}^n|\leq 1$ for $0 < n$ and $1\leq i, j\leq N_x \cdot N_y$.

Proof. First, let us define two infinity norms of $m$-by-$n$ matrix $A$ and $n$-by-$1$ vector $v$. They are respectively defined as follows:

Then, we know that $\|Av\|_{\infty}\leq \|A\|_{\infty}\|v\|_{\infty}$. Together with Lemma 3 and adding at most $2/h^2$ to the operator $\mathcal{L}_n$ as the constraints, it implies that the size of the Laplacian with constraints is bounded by

Next, we assume that $\|\phi^n\|\leq 1$, and discuss the numerical solution's boundedness by mathematical induction. Again from Lemma 3 and the condition $0 < \Delta t \leq h^2/(2d+2)$, we see readily that the linear operator $(I_n-\Delta t \mathcal{L}_c)$ is nonsingular. Note that $\|L_{c}\|_{\infty}\leq 1$, and the spectral radius of $L_{c}$ follows $\rho(L_{c})\leq 1$ since $(\mathcal{L}_{c})_{i, i}+\sum_{j\ne i}|(\mathcal{L}_{c})_{i, j}|\leq 1$ for all $i, j\in[1, N_x \cdot N_y]$. Observe that the infinity norm of $(I_n-\Delta t \mathcal{L}_c)^{-1}$ follows from

Together with $\|\phi^*\|_{\infty}\leq 1$, we have that $|\phi^*_{i, j, k}|\leq \sqrt{ (\phi^*_{i, j, k})^2 +e^{-\frac{2\Delta t}{\epsilon^2}} (1-(\phi^*_{i, j, k})^2)}$ since $(1-(\phi^*_{i, j, k})^2)$ and $e^{-\frac{2\Delta t}{\epsilon^2}}$ are always positive. Therefore, the numerical solution satisfies

□

4.   Numerical tests for the proposed method

Now, our aim is to construct a surface with boundary $\Gamma$ and mean curvature zero at all points. We present the two classical examples and one example containing the certain frame boundary as the benchmark cases. In the case of Scherk's surface, the right angle frames on the top of the surface are the main interest of the present work. In the catenoid as the next case example, we consider the round frames in the top and bottom of the surface.

4.1.   Example $1$: Scherk's surface

First, let us consider the numerical construction of Scherk's surface. We set $128\times 128\times 128$ uniform grid points for the computational domain $\Omega_h: = \{(x_i, y_j, z_k)\in (0, 1)\times(0, 1)\times (0, 1) \leq i, j, k\leq 128\}$. The parameters used are the grid size $N = N_x = N_y = N_z = 256$, spatial step size $h = 1/N$, $\Delta t = h^2/16$, and $\epsilon = 4h/(2\sqrt{2}\tan^{-1}(0.9))$.

The boundary condition as zero constraints is created to form a box frame missing two parallel edges on the top and two parallel edges on the bottom. With the index set $S = \{1, \ldots, 8\}$ and width parameter $a_1 = 20$, we define the constraint domains $\Gamma = \cup_{i\in S} \Gamma_i$ as below:

Then, the zero constraints on $\Gamma$ are given by

For numerical purposes, let us vectorize $\phi_{i, j, k}$ by concatenating the elements vertically, i.e., $\phi_{i, j, k} = \phi_{i+(j-1)N_x+(k-1)N_x N_y}$. We define $R_i: = \sum_{j\ne i}|C_{i, j}|$ for $1\leq i, j\leq N_x N_y N_z$. Then, the constraint matrix $C$ from (4.1) follows $|C_{i, i}|\geq R_i\quad \text{for }1\leq i\leq N_x N_y N_z$, where $C_{i, j}$ denotes the entry in the $i$-th row and $j$-th column. This implies that $C$ is positive semidefinite by the Gershgorin circle theorem.

Figure 3 shows the zero level contours of numerical solutions. Each surface in Figure 3(a), (c) and (d) is represented by the zero level sets. The shaded regions in Figure 3(b) are imposed with $\phi = 0$ implicitly as in the equations (4.1).

The frame boundaries are represented by black dotted lines in Figure 4. Surfaces are also illustrated by the zero level set of solutions. On the right side of each figure, the curved surface is attached due to the geometric constraints. The dotted lines are construction lines, indicating right angles and straight lines. The yellow colored surface is held up by the frame all the time. To numerically validate the numerical security of the frame boundary, we take a closeup look at the right-top corner of the surface in Figure 4: (a) closeup at time $t = 0$, (b) closeup at time $t = 1.831\times 10^{-1}$ and (c) closeup at time $t = 9.766\times 10^{-1}$.

Scherk's surface is one of the plateau's surfaces. It forms the shape of a soap film having the boundary of a square which is bent upward on two opposing sides and downward on the other two sides as in Figure 3(a). This surface is represented by $\{ (x, y, z)\in \mathbb R^3|z = \ln(\cos(y)) -\ln(\cos(x))\}$. A piece of Scherk's surface defined on $-1\leq x \leq 1$ and $-1\leq y \leq 1$ is plotted in Figure 5(a). Comparison of Figure 5(a) and (b) confirms numerically that the proposed scheme seems to provide almost the same result by the analytic formulation. For comparison purposes, we excluded the frame boundary conditions of the curved surface, and plot the surface in Figure 5(b).

4.2.   Example $2$: Catenoid

The catenoid is one of the minimal surfaces, and it is readily formed by a soap film stretched across two wire discs, the planes of which are perpendicular to the line joining their centers. The catenoid is a member of the one-parameter family of surfaces of revolution of the catenary $y = a\cosh(x/b)$, where $a$ and $b$ are constants. The parametric equations for the catenoid are given by

We derived the analytic formulation by the Weierstrass representation in the previous section, but we also note that the catenoid can be generated by the calculus of variations. We give an exposition for the completeness of our demonstration.

Let us begin with the catenary. By assuming a continuous solution in $C^2$, we derive a curve which is called a catenary. We find the minimum surface area of revolution among all the curves joining two points $(x_1, u_1)$ and $(x_2, u_2)$. We can obtain the catenary by the calculus of variations. The surface area of revolution is obtained by

Here, we find $y^{*}(x)\in C^2([x_1, x_2])$ as a minimizer of above the function $\text{Area}(x)$. Let $F(y, \dot{y}) = y\sqrt{1+\dot{y}^2}$, and we know that

Then, the Euler-Lagrange equation implies that

This gives us the first integral:

Then, the equation reduces to

With the arbitrary constants, this differential equation can be solved by

Now, let us consider the numerical construction of the catenoid. We perform this test for numerical comparisons between analytic and numerical surfaces. Solutions from the numerical test are radially symmetric, but we are concerned about the connection with zero interface $\mathcal{I}$ of the solution and the rounded frame boundaries.

On the computational unit cubic domain $\Omega_h = (0, 1)\times (0, 1)\times (0, 1)$, the parameters used are the grid size $N = N_x = N_y = N_z = 256$, spatial step size $h = 1/N$, $\Delta t = h^2/16$, and $\epsilon = 4h/(2\sqrt{2}\tan^{-1}(0.9))$, and the initial data as shown in Figure 6(a) is given by

As the frame boundaries, there are two circular bands positioned on top of each other. The bands are upright, and the faces of the bands face the front. We display the specific frame boundaries in Figure 6(b). The shaded regions are imposed with $\phi = 0$ implicitly in the same fashion of Scherk's surface. We always attach the solution to the circular bands, using them as geometric constraints. We show the numerical solution with 40000 iterations in 6(c). The surface is bent in the middle part of the cylinder. Then, the bending motion lasts no longer in Figure 6(d). We simulate, resulting in the numerical solution after $1.0\times 10^{6}$ iterations at time $t = 0.9537$ when solution converges.

Figure 7 presents the final shape of the zero level-set of the numerical solution. Since the frame boundaries are not part of the minimal surface, we get rid of the frames for schematic definiteness. The Surface from the numerical solution is overlapped by the gray colored analytic surface. Three different (tilted, side, top) perspectives are illustrated in Figure 7. The left portion of the merged surface is constructed using numerical solutions, with the top and bottom being blocked by representation of the zero level-set. As we can see that, all results are consistent with the analytic surface.

4.3.   Example $3$: Two right-angled bands

So far, we have dealt with surfaces with well-known analytical solutions, i.e., Scherk's surface and catenoid. In this example, we construct a surface by using the proposed numerical method for the case where the analytic solution is unknown. Apart from the surfaces that exhibit the periodicity in three directions, i.e., triply periodic minimal surfaces, there are only a few minimal surfaces with ends in differential geometry. Even if such minimal surfaces are known, the majority of them are highly unstable to construct numerically. For this reason, stability is necessary to be stable in order to find minimal surfaces numerically. In this case, let us examine a surface with stability similar to the catenoid but with the ends forming right angles with straight lines.

We now consider two parallel band constraints. Two bands are formed by straight lines at right angles. This is a problem that does not have a known solution. However, in this problem, we can approximate the analytic solution as a numerical solution. There are two bands, one on top of the other. The bands are standing upright, and the surfaces of the bands face the front. This is depicted in Figure 8(a). The grayed parts corresponds to the frame boundaries. We restrict the solution based on these bands as the frame boundary conditions, and observe the evolution of the rest of the surface. The initial datum as shown in Figure 8(b) is given by

The boundary condition as zero constraints is created to form two parallel and angled bands. The frame boundary conditions are divided into two parts, one below, $\Gamma_1, \; \Gamma_2, \; \Gamma_3, \Gamma_4$, and one above, $\Gamma_5, \; \Gamma_6, \; \Gamma_7, \; \Gamma_8$. With the index set $S = \{1, \ldots, 8\}$ and width parameter $a_1 = 20, \; a_2 = 60$, we define the constraint domains $\Gamma = \cup_{i\in S} \Gamma_i$ as below:

Then, the zero constraints on $\Gamma$ are given by

The AC equation is derived from the Ginzburg-Landau (GL) energy functional

The critical points of the above energy functional are of interest for area-minimizing surfaces. By decreasing the GL energy functional (4.3), we can locally minimize the area of each point on the surface. In this regard, there have been several studies on the connection between the Allen-Cahn equation and the minimal surface (^[11,12,25] and references therein). In this third example, we can numerically confirm that the surface obtained by the proposed method is a minimal surface, as it decreases the GL energy functional (4.3).

We discretize the GL energy functional (4.3), and calculate the energy in order to show numerically that the Ginzburg-Landau energy functional decreases. The GL functional can be discretized at time $t^k$ as

Figure 9 shows that the GL energy functional decreases over time $t$, and eventually reaches a point where the decrease stops. The energy is initially high, as illustrated at $t = 0$ in Figure 8(b) and Figure 9; but as it decreases rapidly, the shape of the surface also undergoes a rapid change, as seen in Figure 8(c). From a certain point on, the surface hardly evolves, as seen in Figure 8(d), and the energy also stops decreasing.

5.   Conclusions

The connection of the Allen–Cahn (AC) equation and minimal surfaces has been studied over the past decades, and several numerical schemes have been proposed and developed. However, there is an obstacle to dealing with the fixed frame boundary. Often, the frame boundary makes the numerical scheme unstable and unable to construct minimal surfaces. As far as we know, there is no research work on the frame boundary conditions.

Therefore, we propose a numerical scheme, which is stable and secure for the frame boundaries. We have proven that the AC's Laplacian operator with the frame boundary is still nonsingular, and solutions are bounded by our proposed numerical scheme. In addition, we have demonstrated two classical examples, Scherk's surface and catenoid and one surface for which the analytical solution is unknown. In these cases, we observed the frame boundaries containing right angle and straight lines or rounded frames. As seen before, we have checked that the proposed scheme works well in various frame conditions.

Last, we notice that the construction of some minimal surfaces can be impossible or unstable due to the geometric properties that the surface inherently has. These instabilities are not related to the frame boundaries. However, there is no consensus on the stability of minimal surfaces. We hope that our ongoing research will guide us to a better understanding of the AC equation and minimal surfaces.

Use of AI tools declaration

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

Acknowledgment

This work was supported by Incheon National University Research Grant in 2021-0357.

Conflict of interest

The author confirms that this article content has no conflict of interest.