Statistical analysis of numerical solutions to constrained phase separation problems

1.   Introduction

We consider a phase separation problem that involves the minimization of the functional

where $\epsilon > 0$ and $\beta > 0$ are constant parameters. The integral in Eq (1.1) is taken over a compact surface $M\subset\mathbb{R}^3$, $dS$ is surface area measure, and $\nabla_S$ is a surface gradient operator. The functional $F_p(\phi)$ is a Landau-type free energy that is regularly considered in phase separation problems. Indeed, minimization problems using $F_p(\phi)$ are studied in ^[1,2,3] with the linear conservation constraint

where $\omega \in (0, 1)$ is a conservation constraint parameter and $|M|$ is the surface area of $M$. The problem of minimizing Eq (1.1) subject to Eq (1.2) with $\phi$ in an appropriate space of admissible functions is sometimes called the diffuse interface problem. Using $\Gamma$-convergence methods, one can study a related functional that is defined on certain subsets $\mathcal{E}\subset M$ with $|\mathcal{E}| = \omega |M|$. This latter approach gives rise to a problem known as the sharp interface problem (see, e.g., ^[4] or ^[5]).

The literature concerning relationships between the diffuse interface problem and the sharp interface problem is vast. As in ^[6], we study how numerical solutions arising in the diffuse interface approach compare to those from the sharp interface problem. Let $\chi_\mathcal{E}$ be the indicator function for $\mathcal{E}\subset M$, i.e., $\chi(x, y, z) = 1$ if $(x, y, z)\in M$ and $\chi(x, y, z) = 0$ if $(x, y, z)\notin M$. If $\phi = \chi_\mathcal{E}$, then Eq (1.2) reduces to the area constraint $|\mathcal{E}| = \omega |M|$. In this work, we will compare results obtained from minimizing $F_p$ subject to the linear constraint Eq (1.2) with results obtained from minimizing $F_p$ subject to a nonlinear constraint in the form of

for an appropriate nonlinear function $f(\phi)$. Our particular choice is $f(\phi) = (\phi^2 - 2\phi)^2$, and we note that Eq (1.3) also reduces to $|\mathcal{E}| = \omega|M|$ if $\phi = \chi_\mathcal{E}$. Our motivation for considering such a nonlinear constraint stems from ^[7] wherein a similar function of $\phi$ is utilized in an effort to better capture the sharp interfacial structure of the boundary between two distinct regions of interest. Our intent in using Eq (1.3) then is to "encourage" the formation of a sharper interface between $\mathcal{E}$ and $M\setminus \mathcal{E}$ than what is found when the minimization problem is subject to the linear constraint.

Specific details will be given in Section 2, but let us list our two problems for suitably defined spaces of admissible functions $\mathcal{A}_{LC}$ and $\mathcal{A}_{NLC}$. When the linear constraint is applied, we have

when the nonlinear constraint is applied, we have

In both problems, the conservation constraint is built into the space of admissible functions. Let the Sobolev space $W^{1, 2}(M)$ be the completion of $C_{pw}^1(M)$, i.e., the space of piecewise continuously differentiable functions on $M$, with respect to the norm

As in other work (see, e.g., ^[8]), we define the admissible function spaces

We consider a two-phase system in this work, where the surface $M$ might be thought of as an approximation of a static lipid bilayer membrane. Physical realizations of such systems have been studied by other authors (see, e.g., ^[9]). If distinct phases correspond to different concentrations of two different types of lipids in the membrane, then for certain values of the parameters we expect the lipid that appears in less quantity to coalesce into a patch. It is observed physically (see ^[9] or ^[3]) and predicted analytically (see ^[3]), that such patches will tend to take the shape of a geodesic disk centered near a point on the surface where the Gauss curvature, $K$, is a maximum provided that the linear conservation constraint parameter $\omega$ is small enough. In ^[1], the authors provide numerical explorations on ellipsoids that confirm these observations. Larger patches, corresponding to larger values of $\omega$, are observed to behave like the small ones. In ^[10], the problem is taken up again by considering numerical solutions on Cassinian ovals. When Cassinian ovals are used instead of ellipsoids, one finds that certain large patches deviate from the geodesic disk shape. The authors of ^[10] develop a scheme for determining how closely a patch resembles an actual geodesic disk. We call this scheme the "geodesic protocol", and we describe the approach in Section 4 so that differences in the quality of solutions arising in the linearly constrained problem and the nonlinearly constrained problem can be assessed. The robustness of the aforementioned procedures has clearly been established in previous works. The protocol for determining relevant features of the patch was developed in ^[10], and the general phase separation problem on compact surfaces has been widely studied in both two and three dimensions (see, e.g., ^[1,2,3,5,8]). Our new approach of determining statistically significant differences between the linear and nonlinear constraint methods could be applied to any compact surface previously studied but is restricted to ellipsoids and Cassinian ovals for this work. These methods differ from what is in the phase separation literature because we are completing a full statistical analysis on the difference between the two approaches and not just computing an absolute difference in one measurable quantity after obtaining solutions in each of two different approaches, as in ^[11].

The rest of the paper will proceed as follows. Section 2 contains information about our numerical approach and discretization. Our numerical solutions on ellipsoids and Cassinian ovals are presented in Section 3. Details about our geodesic protocol are given in Section 4. Statistical analysis of our solutions is presented in Section 5, and the paper concludes in Section 6.

2.   Discretization of problem P$_{LC}$ and problem P$_{NLC}$

We use a finite element method for the discretization of Problem P$_{LC}$ and Problem P$_{NLC}$. The finite minimization problems are then processed with MATLAB's fmincon routine. We use fmincon's interior point method with a user-supplied gradient and Hessian. This routine is chosen since it allows for a straightforward implementation of both the linear equality conservation constraint and the nonlinear equality conservation constraint.

The main format of our approach is standard, and it is essentially the same as the approach used in ^[1,2,10]. We opt to use the reference triangle $T = \{(r, s)\; |\; 0\le s\le 1-r, \; 0\le r\le 1\}$ in the $rs$-plane with vertices at $N_1(0, 0), N_2(1, 0),$ and $N_3(0, 1)$. See Figure 1(a).

Let $\Omega\subset\mathbb{R}^2$ be a suitable domain, and let ${\bf{X}}$ be a parameterization of $M$, i.e., ${\bf{X}}:\Omega\rightarrow \mathbb{R}^3$ with ${\bf{X}}(\Omega) = M$. For convenience, suppose $T\subset\Omega$ and set $V_i = {\bf{X}}(N_i)$ for $i = 1, 2, 3$. Then $V_i$ for $i = 1, 2, 3$ are the vertices of a triangular facet $\tau_j\subset M_h$, where $M_h$ is a faceted surface that approximates $M$ with mean grid size $h$. The edges of $\tau_j$ will be denoted by ${\bf{e}}_{12} = V_2-V_1$, ${\bf{e}}_{23} = V_3-V_2$, and ${\bf{e}}_{31} = V_1-V_3$. In this context, it is common to refer to the vertices as nodes. Let $N_h$ be the number of nodes in $M_h$, and let $N_T$ be the number of triangular facets in $M_h$. We consider $M_h = \bigcup_{j = 1}^{N_T} \tau_j$ as an approximation to $M$, and we use Personn's mesh generator (see ^[12]) to produce $M_h$. Figure 1(b) shows $M_h$ for a Cassinian oval with $N_h = 2040$ and $N_T = 4076$.

As in ^[13], the shape functions for linear finite elements are $S_1 = 1-r-s, S_2 = r$, and $S_3 = s$. We approximate a smooth phase function $\phi$ on $M$ with the piecewise linear function

where $\phi_i = \phi(V_i)$ are phase values at the nodes of $M_h$. $\phi_h$ is defined on $M_h$. The degrees of freedom in our minimization problem are the phase values at the $N_h$ nodes in $M_h$. Let $\Phi$ be the vector of the $N_h$ phase values.

To discretize the phase energy, $F_p$, we deal with its two terms separately. For the first term, we approximate

Let $E, F,$ and $G$ be coefficients of the first fundamental form on $M$. From ^[14],

We use the approximations $E\approx {\bf{e}}_{12}\cdot {\bf{e}}_{12}, F\approx {\bf{e}}_{12}\cdot {\bf{e}}_{23},$ and $G\approx {\bf{e}}_{23}\cdot {\bf{e}}_{23}$ in our model. It follows that

Integration over the reference triangle is done exactly to find coefficients of $\phi_1, \phi_2,$ and $\phi_3$ on each triangular facet. The $N_h\times N_h$ matrix $A$ is formed by filling $3\times 3$ submatrices with such coefficients. For the second term in $F_p$, we compute

As above, the integrals are computed exactly. This results in a nonlinear function of the phase values, $N(\Phi)$. Since we are not pursuing a complete parametric study here, we set $\beta = 0.25$ in our code. This choice of $\beta$ is consistent with ^[10]. Setting $F_p^h(\Phi) = \frac{\epsilon^2}{2}\Phi^T A\Phi + N(\Phi)$, we have $F_p(\phi)\approx F_p^h(\Phi)$.

Our discrete variational problem with linear constraint is then

for a $1\times N_h$ linear constraint matrix $A_{eq}$. Specifically, $A_{eq}$ is computed via the approximation

Note that, with appropriate changes arising from our choice of reference triangle, Problem$^h_{LC}$ is the same as ^[2] with $\gamma = 0$ since we are not including an inhibitory term in the present work.

Our discrete variational problem with nonlinear constraint is

for a nonlinear constraint function ceq that codes our nonlinear conservation constraint Eq (1.3). In writing the function, we use $f(\phi) = (\phi^2 - 2\phi)^2$ and the approximation

where the integration over $T$ is carried out explicitly. $ceq(\Phi)$ is then

The surface areas $|M|$ are computed exactly for the ellipsoids. We use the surface area formula in ^[15] as an approximation to $|M|$ for the Cassinian oval surface. See Table 1 for a list of surface areas and other mesh information. Note that values are rounded to four decimal places for convenience. Also, naming schemes for the various geometries are presented in Section 3.

Remark 2.1. We use $f(\phi) = (\phi^2 - 2\phi)^2$ in Eq (2.3). Similar to the approach in ^[7], we choose $f(\phi)$ so that $f(0) = 0$, $f(1) = 1$, and $f'(0) = f'(1) = 0$. It has been our convention that $\phi\approx\chi_\mathcal{E}$ in an attempt to approximate the sharp interface problem. However, we could just as well have $\phi\approx \chi_{M\setminus \mathcal{E}}$. In this case, one might elect to use $f(\phi) = (\phi^2-1)^2$ in Eq (2.3) with $f(0) = 1$ and $f(1) = 0$. In future work, one might consider such other choices for $f$.

While fmincon is able to handle inequality constraints and bounds on $\phi$, i.e., constants $lb$ and $ub$ such that $lb\le \phi\le ub$, we do not employ those options. Our problem does not contain any inequality constraints, and we find that MATLAB is better able to return "separated" solutions when we set $lb = ub = []$. After MATLAB completes its minimization, numerical solutions are checked against Def. 3.1 to ensure that appropriate upper and lower bounds are satisfied.

3.   Numerical results

We provide an overview of our numerical solutions in this section. We refer to our solutions as patches. We only wish to consider patches that are suitably well-separated. We adopt a definition from ^[10].

Definition 3.1. (a) A strongly separated solution to Problem P$^h_{LC}$ is a collection of values $\phi_h^* = \; \{ \phi_{h, 1}^*, \phi_{h, 2}^*, \ldots, \phi_{h, N_h}^* \}$ such that $F^h_p(\phi_h^*)\le F^h_p(\Phi)$ for all $\Phi\in\mathbb{R}^{N_h}$ with $A_{eq}\Phi = \omega |M|$ and

(b) A strongly separated solution to Problem P$^h_{NLC}$ is a collection of values $\phi_h^* = \{ \phi_{h, 1}^*, \phi_{h, 2}^*, \ldots, \phi_{h, N_h}^* \}$ such that $F^h_p(\phi_h^*)\le F^h_p(\Phi)$ for all $\Phi\in\mathbb{R}^{N_h}$ with $ceq(\Phi) = 0$ and

Note that our definition of a strongly separated solution is essentially the definition of a "$\delta_0$ strongly separated solution" from ^[10] with $\delta_0 = 0.1$ and the additional restrictions that $\max\{\phi_h^*\}\le 1.02$ and $\min\{\phi_h^*\}\ge -0.02$. These additional restrictions are used in lieu of a priori bounds $lb$ and $ub$. The restrictions ensure that strongly separated solutions do not contain maximum and minimum phase values that deviate too far from 0 and 1. We can now be more specific about what is meant by a "small" patch and a "large" patch. A small patch is a strongly separated solution with $0 < \omega\le 0.25$, and a large patch is a strongly separated solution with $0.25 < \omega\le 0.5$.

For some choices of the parameters, fmincon is unable to compute strongly separated solutions. It will sometimes happen that fmincon returns a separated solution with phase values larger than 1.02 and/or smaller than $-0.02$. We do not use such solutions. Other times, MATLAB terminates with a constant solution $\phi_h^* \approx \phi_{const}$. In this case, when the linear conservation constraint is imposed, Eq (1.2) implies that $\phi_{const} = \omega$. It then follows that $F_p(\phi_{const}) = \beta \omega^2(1-\omega)^2|M|$. When the nonlinear conservation constraint is imposed in such situations, Eq (1.3) implies that $\phi_{const} = G(\omega)$ for a suitable constant function $G$ that solves

It then follows that $F_p(\phi_{const}) = \beta G(\omega)^2(1-G(\omega))^2|M|$. As an example, when $\omega = 0.25$, $\epsilon = 0.5$, and $\beta = 0.25$, we find that $\phi_{const}\approx 0.2929$ and $F_p(\phi_{const})\approx 0.3720$ for our oblate ellipsoid with $|M|\approx 34.6875$. This is in good agreement with the MATLAB returned values of $\phi_{const}\approx 0.2930$ and $F_p\approx 0.3718$. Henceforth, we focus solely on strongly separated solutions.

We seek to answer the question "Does using the nonlinear constraint yield better results?" Better results can be interpreted in many different ways. We pursue two directions to this end. We will see that using the nonlinear constraint allows MATLAB to compute strongly separated solutions for parameter values that do not yield strongly separated solutions with the linear constraint. The nonlinear constraint is better in this case in terms of providing a means to compute more strongly separated solutions. Secondly, we will perform statistical hypothesis tests in Section 5 in order to assess if there are statistical differences in the strongly separated solutions that are computed with the nonlinear constraint compared to the solutions that are computed with the linear constraint.

In solving the constrained minimization problems, we use the default values for the Lagrange multipliers that are provided by MATLAB. While it is possible to provide physical interpretations of Lagrange multipliers in a given context when $F_p$ is viewed appropriately, we do not explore this direction in the current work. Such explanation would deviate too far from our main objective of performing statistical analysis of solution quality, regardless of a particular physical meaning. Preliminary observations of Lagrange multiplier values from the linearly constrained and nonlinearly constrained cases seem to suggest little difference between the two, but we leave a more careful analysis for a future study. Additional future work considering physical interpretations of the Lagrange multiplier arising from the nonlinear constraint in comparison to known physical interpretations of the Lagrange multiplier arising from the linear constraint could also be pursued.

In the remainder of this section, we present strongly separated solutions on ellipsoids and Cassinian ovals. In order to compute a strongly separated solution, we input an initial phase configuration into MATLAB of the form

Note that we use the same naming convention, with appropriate corrections in notation, when the piecewise function is defined by restricting $x$ or $y$. Figure 2 shows the initial phase configuration $\phi_{h, x > 1}^0$ that is used for our Cassinian oval solutions.

3.1.   Summary of results for ellipsoids

We solved the phase separation problems on ellipsoids

We will identify each ellipsoid with the naming scheme $(a, b, c)$ ellipsoid. Our results are from a $(1, 1, 2)$ prolate ellipsoid and a $(2, 2, 1)$ oblate ellipsoid. All of the solutions on the $(1, 1, 2)$ prolate ellipsoid are computed with an initial phase configuration $\phi_{h, z > 1.75}^0$. All of the solutions on the $(2, 2, 1)$ oblate ellipsoid are computed with an initial phase configuration $\phi_{h, y > 1.75}^0$.

Solutions are denoted by $(\omega, \epsilon)$ pairs. For example, we found a $(0.05, 0.04)$ strongly separated solution on the $(1, 1, 2)$ prolate ellipsoid. This means that we selected the parameters $\omega = 0.05, \epsilon = \; 0.04, a = b = 1, c = 2$, and $\phi_{h, z > 1.75}$ and ran our minimization program. MATLAB returned a phase configuration $\phi_h^*$ satisfying Def. 3.1. In this work, we do not claim to have conducted an exhaustive search of the $\omega\epsilon$-plane. Instead, the solutions on the ellipsoids were found with a search over $0.05\le \omega\le 0.5$ with step size 0.05 and $0.04\le \epsilon\le 0.12$ with step size 0.01 in the small $\omega$ case and $0.05\le \epsilon\le 0.3$ with step size 0.05 in the large $\omega$ case. In Table 2, example solutions are given for a $(1, 1, 2)$ prolate ellipsoid. "No Change" means that the same set of $\epsilon$ values that worked for the given $\omega$ in the linearly constrained problem worked for the nonlinearly constrained problem.

Of note, when $\omega = 0.05$, we obtain four additional strongly separated solutions when the nonlinear conservation constraint is used. When $\omega = 0.1$, we obtain three additional strongly separated solutions when the nonlinear conservation constraint is used. The $(0.15, 0.04)$ and $(0.2, 0.04)$ solutions to Problem P$^h_{NLC}$ had $\min\{\phi_h^*\} < -0.02$. Oddly, the $(0.5, 0.2)$ solution to Problem P$^h_{NLC}$ resulted in $\phi_h^*\approx\phi_{const}\approx1.5$.

In Table 3, example solutions are given for a $(2, 2, 1)$ oblate ellipsoid. These solutions were found using the same search in the $\omega\epsilon$-plane as was used for the presented prolate ellipsoid solutions. Additionally, we computed the $(0.3, 0.04)$ solution in order to have a sample size of at least thirty solutions for our statistical analysis in the large solution regime. Of note, when $\omega = 0.05$, we obtain five additional strongly separated solutions when the nonlinear conservation constraint is used. When $\omega = 0.1$, we obtain one additional strongly separated solution when the nonlinear conservation constraint is used.

3.2.   Summary of results for Cassinian ovals

We solved the phase separation problems on a Cassinian oval

with shape parameters $a = 1$ and $c = \sqrt{\frac{1+\sqrt{17}}{4}}$. As in ^[10], these constants are chosen so that the Gauss curvature at $x = y = 0$ is equal to the Gauss curvature at $z = 0$. The Cassinian oval geometry is interesting to consider since it is qualitatively similar to a discocyte Willmore surface. Additionally, it is an example of a geometry that contains regions with $K > 0, K = 0$, and $K < 0$, unlike the ellipsoids that have served as benchmark geometries in previous work. Although the Cassinian oval geometry contains three distinct regions of Gauss curvature, each strongly separated solution is centered at a point in the positive curvature region. All of the solutions on the $\left(1, \sqrt{\frac{1+\sqrt{17}}{4}} \right)$ Cassinian oval are computed with an initial phase configuration $\phi_{h, x > 1}^0$.

In Table 4, example solutions are given for a $\left(1, \sqrt{\frac{1+\sqrt{17}}{4}}\right)$ Cassinian oval. See Figure 3 for an image of the $(0.5, 0.15)$ strongly separated solution to Problem P$^h_{NLC}$. We note that the Table 4 solutions were found with a search considering only $0.05 \leq \omega \leq 0.5$ with a step size of $0.05$ and $0.02 \leq \epsilon \leq 0.1$ with a step size of $0.02$. With $\epsilon\in\{0.11, 0.12, 0.13, 0.15\}$, the $(0.45, \epsilon)$ and $(0.5, \epsilon)$ solutions to Problem P$^h_{NLC}$ were computed so that the total sample size of such large solutions would be at least thirty. Of note, when $\omega = 0.05$, we obtain five additional strongly separated solutions when the nonlinear constraint is used. We obtain two additional strongly separated solutions when $\omega = 0.1$ and when $\omega = 0.15$. The $(0.4, 0.05)$, $(0.45, \epsilon)$, and $(0.5, \epsilon)$ for $\epsilon\in\{0.04, 0.05, 0.06, 0.07\}$ parameter pairs failed to produce strongly separated solutions to Problem P$^h_{NLC}$.

In summary, despite not conducting an exhaustive search of the $\omega\; \epsilon$ -plane for strongly separated solutions, we see that the use of the nonlinear conservation constraint allows for strongly separated solutions to be computed more readily for the smallest $\omega$ values. This is an important benefit to using the nonlinear constraint since small patches are typically more difficult to compute than larger patches. Generally speaking, the required computing time is larger when the nonlinear constraint is applied. As such, using the linear constraint might be preferred for larger patches. Furthermore, our computed large solutions suggest that the solution space to Problem P$^h_{NLC}$ might be smaller than the solution space to Problem P$^h_{LC}$ in the large regime.

4.   Geodesic protocol

Once a solution to Problem P$^h_{LC}$ or Problem P$^h_{NLC}$ has been computed, the solution is checked against Def. 3.1. to determine whether or not it is a strongly separated solution. Solutions that are deemed strongly separated then enter into the geodesic protocol that we describe in this section. The center of the patch $(x_c, y_c, z_c)$ is determined first, based on the particular surface geometry. Let $(\overline{X}, \overline{Y}, \overline{Z})$ be the patch center of mass coordinates that MATLAB computes by direct integration of the formulae

We set $(x_c, y_c, z_c) = t_c(\overline{X}, \overline{Y}, \overline{Z})$ where $t_c$ is determined for the ellipsoids and Cassinian ovals by plugging $t_c(\overline{X}, \overline{Y}, \overline{Z})$ into Eqs (3.1) and (3.2), respectively.

After determining $(x_c, y_c, z_c)$, all nodes with corresponding $\phi_h^* > 0.75$ are labeled as members of the patch, i.e, the patch is $\mathcal{E} = \{(x, y, z)\in M_h\; |\; \phi_h^*(x, y, z) > 0.75\}$. From here, the radius of the patch, $\rho$, is found by employing a shooting method together with MATLAB's ode45 to approximate geodesic curves on the surface. $\rho$ is defined as the length of the longest geodesic connecting the center of the patch to a node in the patch.

Four statistics of interest can now be computed for each solution. Let $N_\mathcal{E}$ be the number of nodes in $\mathcal{E}$. As in ^[10], we identify nodes that are in violation of a decision criteria. Specifically, such nodes in violation are vertices $(x, y, z)\in V$ where

and $d_g({\bf{X}}; {\bf{Y}})$ is the geodesic distance between points ${\bf{X}}, {\bf{Y}}\in M_h$. In the definition of $V$, the tolerance $0.03$ is chosen arbitrarily, but it suffices for our purposes of identifying nodes that are close to $\partial \mathcal{E}$. Four statistics of interest are listed as follows.

1) The number of nodes in violation, denoted by $N_V$.

2) The mean phase value of the nodes in violation, denoted by $\bar\Phi$.

3) The proportion $N_V/N_\mathcal{E}$.

4) $\delta = \min\left\lbrace |d_g((x, y, z);(x_c, y_c, z_c)) - \rho|\; \big|\; (x, y, z)\in V\right\rbrace$

Patch nodes, nodes in violation, and $\rho$ are illustrated in the schematic Figure 4. These statistics were used in ^[10] in order to measure patch deviation from a geodesic disk shape.

5.   Statistical analysis

With an aim toward quantifying differences between solutions to Problem P$^h_{LC}$ and Problem P$^h_{NLC}$ we will test for a mean difference by considering hypotheses in both the small $\omega$ and large $\omega$ regimes. We restrict our analysis to the oblate ellipsoid and Cassinian oval geometries in order to best demonstrate our methods.

The hypothesis tests are conducted for mean differences in eight quantities: $N_V$, $\bar\Phi$, $P = \frac{N_V}{N_\mathcal{E}}$, $\delta$, $N_\mathcal{E}$, $\rho$, $F_p(\phi_h^*)$, and total computation time $\mathcal{T}$. We implement the Student's $t$-test for our analysis and, in alignment with standard practices, accept $p$-values of 0.05 and less as statistically significant. This $p$-value measures the probability of our calculated value appearing given that the null hypothesized difference of zero is true. The smaller the $p$-value, the more statistically significant the result since it represents a smaller probability of that result occurring by chance under the null hypothesis. Since we are conducting $t$-tests, we will report a sample mean difference and sample standard deviation of the differences for each test. Sample sizes are needed in order to conduct the hypothesis tests. As is standard, randomly selected samples of size at least thirty should be used. While a large enough sample size is used for much of the work that follows, in some instances sample sizes are slightly less than thirty. Small solutions are particularly difficult to obtain when the linear constraint is used. Table 4 shows that our sample size for small Cassinian oval solutions is $n = 24$. Furthermore, we are unable to select random samples since we are limited by the strongly separated solutions that we are able to compute. The sample mean difference and sample standard deviation of the differences are computed by pairing the data. For example, with the oblate ellipsoid geometry, there is a $(0.05, 0.04)$ strongly separated solution to Problem P$^h_{LC}$ and a $(0.05, 0.04)$ strongly separated solution to Problem P$^h_{NLC}$.

A pair of solutions represents one datum in our analysis. In other words, there is one difference in each of our eight quantities of interest for each single pair of solutions. Differences in quantities are formed by subtracting the quantity from the solution to Problem P$^h_{NLC}$ from the quantity from the solution to Problem P$^h_{LC}$. Figure 5 shows a pair of oblate ellipsoid solutions. For this pair, there is a slight but noticeable difference in patch size. In Figure 5(a), the solution has $0.0139\le \phi_h^*\le 1.0143$ and $F_p(\phi_h^*) = 0.0225$. In Figure 5(b), the solution has $-0.0052\le \phi_h^*\le 1.0011$ and $F_p(\phi_h^*) = 0.0243$. Figure 6 shows a pair of Cassinian oval solutions. Note that while the two images are quite similar, the solutions are not identical. The $(0.3, 0.04)$ solution to Problem P$^h_{LC}$ has $-0.0015\le \phi_h^*\le 1.0065$ while the solution to Problem P$^h_{NLC}$ has $-0.0051\le \phi_h^*\le 1.003$. The free energies of these solutions differs only slightly. The solution to Problem P$^h_{LC}$ has $F_p = 0.0390$ while the solution to Problem P$^h_{NLC}$ has $F_p = 0.0389$.

We find both statistically significant and insignificant differences in our solutions. In this section, we highlight the statistically significant differences. The hypothesis tests yielding insignificant results are included in the Appendix for completeness. At the standard 5% significance level, there is a significant difference in mean $N_\mathcal{E}$ and mean total computation time for both geometries in both size regimes. For small and large oblate ellipsoid solutions and small Cassinian oval solutions, we find significant results at the 5% significance level for mean differences in $\rho, \frac{N_V}{N_\mathcal{E}},$ and $F_p(\phi_h^*)$. The notation $\mu_{*, *}$ will be used to represent a population mean difference. Our notational convention for the subscript is quantity, regime size. For example, the subscript $N_\mathcal{E}, s$ pertains to the number of patch nodes for solutions from the small $\omega$ region $0 < \omega\le 0.25$. On the other hand, the subscript $N_\mathcal{E}, l$ signifies that the quantity refers to solutions from the large $\omega$ region $0.25 < \omega\le 0.5$.

While not a feature of a solution patch, total program run time might be of interest to those researchers with limited computational resources. For example, nearly all of the ellipsoid solutions in this work were computed on a single core personal laptop. The computer resources were taxed, and it was not possible to compute solutions on finer meshes that one could consider on a more powerful machine. The first two hypotheses are for a mean difference in the total program run time:

The total program run time is the sum of the run time for the minimization of $F_p$ and the run time for the geodesic protocol. Differences in total run time are computed for pairs of solutions produced by a single machine. As such, the $\omega = 0.5$ oblate ellipsoid solutions are excluded since they were found on a different machine from the other oblate ellipsoid solutions. Table 5 shows relevant statistics for $HI_s$ and $HI_l$ on our geometries. The reported sample mean and sample standard deviation values are in seconds. Our $p$-values and negative sample mean differences demonstrate that solutions computed with the nonlinear constraint take significantly longer to run than solutions computed with the linear constraint. For those with limited computational power, our statistical results should be considered since we find insignificant differences in solution patches for some regime sizes and geometries. Solving Problem P$^h_{LC}$ is significantly faster without a significantly different solution in such situations.

The next two hypotheses are for a mean difference in $N_\mathcal{E}$:

Table 6 shows relevant statistics for $HII_s$ and $HII_l$ on our geometries. For both geometries, small patches contain more nodes when the nonlinear constraint is used. Figure 5 clearly demonstrates this significant effect. By contrast, our results show that, for both geometries, large patches contain more nodes when the linear constraint is applied.

The third set of hypotheses are for a mean difference in $\rho$:

Table 7 shows relevant statistics for $HIII_s$ and $HIII_l$ on our geometries. The large regime solutions on the Cassinian oval only yield a significant difference at the 10% significance level. On both geometries, small patches have a smaller geodesic radius when the linear constraint is used. By contrast, large patches on both geometries have a smaller geodesic radius when the nonlinear constraint is applied. This is consistent with the results from $HII_k$.

Remark 5.1. A significant mean difference in $\rho$ is particularly important to note for small solutions. As an example, the geodesic radius of our $(0.05, 0.04)$ oblate ellipsoid solution to Problem P$^h_{LC}$ is $0.5990$. The geodesic radius of our $(0.05, 0.04)$ oblate ellipsoid solution to Problem P$^h_{NLC}$ is $0.7563$. Using $\pi\rho^2 \approx \omega|M|$, we compute $0.05|M|\approx 1.7344$. In the nonlinear case, we have $|\mathcal{E}| \approx 1.7970$. In the linear case we have $|\mathcal{E}|\approx 1.1272$. These values correspond to percentage errors of 3.61% and 35.01% in the nonlinear and linear cases, respectively. Using the nonlinear constraint in this instance provides a solution that is much closer to the size of a solution to the sharp interface problem.

The next two hypotheses are for a mean difference in $F_p$:

Table 8 shows relevant statistics for $HIV_s$ and $HIV_l$ on our geometries. The large regime solutions on the Cassinian oval do not yield a significant mean difference in free energy. This can be seen in our discussion of the (0.3, 0.04) solution shown in Figure 6. For both size regimes on the oblate ellipsoid and the small size regime on the Cassinian oval, we see that the free energy is significantly larger when the nonlinear constraint is used. That said, the difference could be slight. For example, the mean difference only accounts for approximately 3.6% of $F_p(\phi_h^*)$ for the $(0.05, 0.04)$ solution to Problem P$^h_{LC}$ on the oblate ellipsoid.

The last two hypotheses in this section are for a mean difference in $P = \frac{N_V}{N_\mathcal{E}}$:

Table 9 shows relevant statistics for $HV_s$ and $HV_l$ on our geometries. We notice an apparent geometry effect upon inspection of the $p$-values. Solutions on the oblate ellipsoid show significant differences in the proportion $\frac{N_V}{N_\mathcal{E}}$. This can be viewed as a consequence of our result that small solutions tend to have fewer patch nodes when the linear constraint is applied compared to when the nonlinear constraint is used. Table A1 in the Appendix shows that there are insignificant mean differences in $N_V$ for all geometries and all size regimes. If $N_V$ is neglected, smaller $N_\mathcal{E}$ values produce larger proportions $\frac{N_V}{N_\mathcal{E}}$. Small solutions on the Cassinian oval only show a significant mean difference at the 10% significance level.

Remark 5.2. We see from Table A1 that there is not a statistically significant mean difference in the number of nodes in violation for solutions to Problem $P^h_{LC}$ and solutions to Problem $P^h_{NLC}$. If a refined patch boundary is interpreted to mean that there are fewer nodes in violation, then we conclude that there is not evidence that solutions to one problem have a more refined patch boundary than solutions to the other problem.

Remark 5.3. Upon closer inspection of our data, one can find more differences between solutions for the smallest $\omega$ values. Some such differences will be presented in the Appendix. As an example, for the smallest ten Cassinian oval solution pairs, there is a significant mean difference in $\frac{N_V}{N_\mathcal{E}}$. The $p$-value of $0.044$ is computed via a randomization method in StatKey (see ^[16]) due to the small sample size of $n = 10$. This result and others like it lead us to suspect that a further refinement of our regime sizes, perhaps based on statistical differences in solutions, might yield interesting results in a future work.

6.   Conclusions

We present a thorough investigation into benefits of using a nonlinear constraint instead of a linear constraint in a well-known phase separation problem. We find that the nonlinear constraint seems to enlarge the solution space in the sense that more strongly separated solutions can be obtained on a fixed mesh size for small $\omega$ values. We identify a small and large $\omega$ regime, and we compute many strongly separated solutions within each regime. Such work is done for two different geometries: ellipsoids and Cassinian ovals. Our subsequent statistical analysis of the solutions shows that there are significant statistical differences between solutions to the linearly constrained problem and the nonlinearly constrained problem. We employ a geodesic protocol to generate statistics that can be used to quantify mean differences between strongly separated solutions to Problem P$^h_{LC}$ and Problem P$^h_{NLC}$. Statistically significant mean differences in total computation time and quantities related to patch size are highlighted for oblate ellipsoids and Cassinian ovals.

Acknowledgments

Converse University provided a Summer Research Grant stipend through the Faculty Development Committee for Amanda Mangum's work toward this study during Summer 2021.

Conflict of interest

The authors declare there is no conflicts of interest.

Appendix

Additional Hypothesis Tests

Three sets of hypothesis tests yield insignificant results that we present here for completeness. Additionally, some hypothesis tests only demonstrate significant results at the 10% significance level, for a one-sided test, or for certain subsets of our data. These results are also presented in this Appendix.

We conducted hypothesis tests for a mean difference in number of nodes in violation:

Table A1 contains relevant statistics for $HVI_s$ and $HVI_l$ on our geometries. As noted in Section 5, Table A1 shows that there is not a significant mean difference in $N_V$ for solutions to the linearly constrained problem and the nonlinearly constrained problem. This suggests that, on average, either $N_V$ is not able to distinguish a sharper interface in the nonlinearly constrained case or that solution patches from Problem P$^h_{NLC}$ do not have a sharper boundary than those from Problem P$^h_{LC}$ for a fixed mesh.

Remark A.1. If we only consider the first ten paired data values in the small regime on the oblate ellipsoid, we compute $p$-value $\approx 0$. In this case, we have significant data that, on average, the number of nodes in violation for solutions to Problem P$^h_{NLC}$ is 4.5 larger than the number of nodes in violation for solutions to Problem P$^h_{LC}$. This $p$-value is computed with a randomization method in StatKey (see ^[16]) due to the small sample size $n = 10$.

The next two hypotheses are for a mean difference in the mean phase value of nodes in violation:

Table A2 shows relevant statistics for $HVII_s$ and $HVII_l$ on our geometries. While the two-sided test for a mean difference is not significant for the small oblate ellipsoid solutions, we do see significant results at a 10% significance level in the one-sided case ($p$-value $\approx 0.082$). This lends some support to the nonlinearly constrained problem being better since it appears that solutions to Problem P$^h_{LC}$ exhibit larger $\bar\Phi$, on average.

Remark A.2. If one considers only the first ten paired values in the small regime for the oblate ellipsoid, the significant mean difference in $\bar\Phi$ is even more pronounced since the two-sided test in this case yields a $p$-value $\approx 0.084$.

The last two hypotheses are for a mean difference in $\delta$:

Table A3 shows relevant statistics for $HVIII_s$ and $HVIII_l$ on our geometries. There is little to no evidence of a mean difference in $\delta$.

Nomenclature

We provide a summary of our nomenclature for the convenience of the reader. Tables A4 and A5 contains notation first appearing in Sections 1 and 2.

Tables A6 and A7 contain notation first appearing in Sections 3–5. The subscript $k\in\{s, l\}$ refers to either the small $\omega$ regime ($k = s$) or the large $\omega$ regime ($k = l$).