1. Introduction

A major challenge in many applied mathematical problems, such as those arising from materials science, is the need to compute large systems of equations. When the problems are cast in terms of ordinary differential equations (usually by imposing some symmetry), this is equivalent to studying equations with higher order derivatives. If one could approximate these higher order differential equations by second order equations, it would be a considerable simplification. Many problems in applied mathematics are associated with interfaces or moving boundaries, and surface tension often plays an important role in those problems [11] as a stabilizing agent for interfaces that would otherwise be highly unstable (e.g., dendritic behavior [2,9,4]). In a recent paper, quantum field theoretic renormalization methods were used to approximate higher order ODEs by second order counterparts [6]. A related issue is whether one can approximate solutions of one higher order equation with those of another with different, e.g., smaller coefficients, that can be used in conjunction with the renormalization methods cited above.

In this paper, we consider a prototype free energy of the phase field type that has been studied in several papers including [5,10], and the resulting differential equations that follow from minimization of this free energy. In particular we show that the suface tension σ as a function of rescaled coefficients satisfies Clairaut's equation that can be solved by the method of charactertics. This means that if one knows the surface tension for a suitable surface Γ in the coefficient space, then one can solve for σ for all values of the parameters.

This methodology can be viewed as complementary to those of [6] involving reduction in the order of the differential equation.

---

2. The Free Energy and Higher Order Phase Field Equations.

We consider the simplest free energy functional that has the necessary features (see [10,5]). Let $\phi\in C^{N}\left (\mathbb{R}\right)$ and let the $j^{th}$ derivative of $\phi$ in $x$ by denoted by either $D^{j}\phi$ or $\phi^{\left (j\right) }.$ Let $\phi$ satisfy the boundary conditions

$\phi^{\left( j\right) }\left( \pm\infty\right) = 0$

(2.1)

and define $W$ as the standard double well potential with minima at $\phi : = \pm\phi_{0},$ e.g.,

$W\left( \phi\right) : = \left( \phi^{2}-1\right)$

A free energy, $\mathcal{F}\left[\phi\right]$, is defined by

$\mathcal{F}\left[\phi\right] : = \int_{\mathbb{R}}F\left( x, \phi\left( x\right), ..., D^{N}\phi\left( x\right) \right) dx \\ F\left[x, \phi\left( x\right), ..., \phi^{\left( N\right) }\left( x\right) \right] = \frac{1}{2}\sum_{j = 1}^{N}\left( -1\right) ^{j+1}% c_{2j}\left\{ D^{j}\phi\left( x\right) \right\} ^{2}+W\left( \phi\left( x\right) \right) .$

(2.2)

This free energy arises from an averaging process in the statistical mechanics of a two state system, where $\phi$ represents an order parameter that makes a transition from the lower energy phase (e.g., solid) at $-1$ (more generally $\phi_{-}$) to $+1$ (more generally $\phi_{+}).$

The higher order phase field equations are obtained from this free energy. For a smooth test function that also satisfies the same boundary conditions as $\phi,$ the functional derivative is evaluated as

$\begin{gathered} \frac{\partial }{{\partial \varepsilon }}\mathcal{F}\left[{\phi + \varepsilon \eta } \right]{|_{\varepsilon = 0}} = \frac{\partial }{{\partial \varepsilon }}\% \int_\mathbb{R} F \left[{x, \phi \left( x \right) + \varepsilon \eta \left( x \right), ..., {\phi ^{\left( N \right)}}\left( x \right) + \varepsilon {\eta ^{\left( N \right)}}\left( x \right)} \right]dx{|_{\varepsilon = 0}} \hfill \\ = \int_\mathbb{R} {\sum\limits_{j = 1}^N {{c_{2j}}} } {\left( { - 1} \right)^{j + 1}}{D^j}\% \phi \left( x \right){D^j}\eta \left( x \right) + {W^\prime }\left( {\phi \left( x \right)} \right)\eta \left( x \right)dx. \hfill \\ \end{gathered}$

Upon performing integration by parts and setting this expression to zero, one has

$0 = \frac{\partial}{\partial\varepsilon}\mathcal{F}\left[\phi+\varepsilon \eta\right] |_{\varepsilon = 0} = \left\{ -\sum_{j = 1}^{N}c_{2j}D^{2j}\phi\left( x\right) +W^{\prime}\left( \phi\left( x\right) \right) \right\} \eta\left( x\right) dx.$

Since this is true for all test functions $\eta,$ one has the ODE,

$-\sum_{j = 1}^{N}c_{2j}D^{2j}\phi\left( x\right) +W^{\prime}\left( \phi\left( x\right) \right) = 0.$

(2.3)

The surface tension in physical terms is usually described in (in arbitrary physical dimension) as the difference in a small cylindrical volume along the surface, normalized by the cross-sectional area of the cylinder. The difference between the free energy obtained by integrating over the cylinder minus the average of the free energy of the two phases is normalized by the area. With $\mathcal{F}_{cyl}\left[\phi\right]$ denoting this free energy one can write (see e.g., [3] or [5])

$\sigma\ \sim\frac{\mathcal{F}_{cyl}\left[\phi\right] -\frac{1}{2}\left\{ \mathcal{F}_{cyl}\left[\phi_{+}\right] +\mathcal{F}_{cyl}\left[\phi _{-}\right] \right\} }{Cross\ Sectional\ Area}.$

(2.4)

In the context of our one-dimensional analysis using a symmetric $W,$ we note that the two terms involving the pure phases vanish, and the integral over this cylinder can be regarded as one-dimensional after the division.

Recalling the comment above that the $W\left (\phi_{\pm}\right)$ terms vanish, it is evident that in one-dimension (or a physical setting with this symmetry) the definition can be interpreted as the free energy of the transition layer solution, $\phi.$ Thus we can write the mathematical definition as follows.

Definition. Given a set of non-negative coefficients, let $\phi$ be a solution to (2.3) subject to boundary conditions (2.1), the surface tension $\sigma\left (c_{2}, ..., c_{2N}\right)$ is defined by

$\sigma\left( c_{2}, ..., c_{2N}\right) = \mathcal{F}\left[\phi\right] .$

(2.5)

Proposition. The surface tension defined by (2.5) can be expressed as

$\sigma\left( c_{2}, ..., c_{2N}\right) = \sum_{j = 1}^{N}\left( -1\right) ^{j+1}jc_{2j}\left\Vert D^{j}\phi\right\Vert _{L^{2}\left( \mathbb{R}\right) }^{2}$

(2.6)

where $\left\Vert f\right\Vert _{L^{2}\left (\mathbb{R}\right) }$ is the usual $L^{2}$ norm, i.e., $\left\Vert f\right\Vert _{L^{2}\left (\mathbb{R}\right) }: = \int_{\mathbb{R}}\left\vert f\right\vert ^{2}dx.$

Remarks. Mathematically, the surface tension is well-defined so long as one has a solution to the equation (2.3). The connection between the two definitions is easily understood (see [7,5]) in light of the calculations in the proof below, which was presented in [7].

Proof. First we derive an identity by muliplying (2.3) by $D\phi$ and integrating over $\left (-\infty, x\right)$

$\int_{-\infty}^{x}-\sum_{j = 1}^{N}c_{2j}\left\{ D^{2j}\phi\left( z\right) \right\} \left\{ D\phi\left( z\right) \right\} +D\phi\left( z\right) W^{\prime}\left( \phi\left( z\right) \right) dz = 0.$

(2.7)

The term involving $W^{\prime}$ is an exact differential, $\frac{d} {dz}W\left (\phi\left (z\right) \right).$ Noting the identity [7]

$D^{2j}\phi D\phi = D\left\{ \frac{\left( -1\right) ^{j-1}}{2}\left( D^{j}\phi\right) ^{2}+\left( -1\right) \sum_{k = 1}^{j-1}\left( -1\right) ^{k-1}D^{2j-k}\phi D^{k}\phi\right\}$

one observes that the left hand side is also an exact identity, yielding,

$\sum_{j = 1}^{N}c_{2j}\left\{ \frac{\left( -1\right) ^{j-1}}{2}\left( D^{j}\phi\right) ^{2}+\left( -1\right) \sum_{k = 1}^{j-1}\left( -1\right) ^{k-1}D^{2j-k}\phi D^{k}\phi\right\} = W\left( \phi\right).$

Note that all terms vanish at $-\infty$ due to the boundary conditions.

Now we integrate this expression over $\left (-\infty, \infty\right)$ and obtain after integrating by parts $j-k$ times in the second term on the left hand side

$\begin{gathered} \sum\limits_{j = 1}^N {\left\{ {{c_{2j}}\frac{{{{\left( { - 1} \right)}^{j - 1}}}}{2}{{\left\| {{D^j}\phi } \right\|}^2} + {c_{2j}}\sum\limits_{k = 1}^{j - 1} {{{\left( { - 1} \right)}^{k - 1}}} {{\left( { - 1} \right)}^{j - k}}{{\left\| {{D^j}\phi } \right\|}^2}} \right\}} \hfill \\ = \int_{ - \infty }^\infty W \left( {\phi \left( x \right)} \right)dx. \hfill \\ \end{gathered}$

Simplifying this expression yields

$\sum_{j = 1}^{N}\left( j-\frac{1}{2}\right) \left( -1\right) ^{j-1}% c_{2j}\left\Vert D^{j}\phi\right\Vert ^{2} = \int_{-\infty}^{\infty}W\left( \phi\left( x\right) \right) dx.$

(2.8)

Using this identity, the free energy (2.2) and the expression (2.5), $\sigma$ can be written as

$\begin{gathered} \sigma \left( {{c_2}, ..., {c_{2N}}} \right)\; = \;\int_{ - \infty }^\infty {\frac{1}{2}} \sum\limits_{j = 1}^N {{{\left( { - 1} \right)}^{j + 1}}} {c_{2j}}{\left\{ {{D^j}\phi \left( x \right)} \right\}^2} + W\left( {\phi \left( x \right)} \right)dx \hfill \\ = \sum\limits_{j = 1}^N {{{\left( { - 1} \right)}^{j + 1}}} j{c_{2j}}{\left\| {{D^j}\% \phi } \right\|^2}, \hfill \\ \end{gathered}$

(2.9)

where the last expression is obtained by substituting (2.8) for the $W$ term.///

---

3. Surface tension and dimensional analysis.

Within this mathematical setting we have been using reduced dimensional parameters. The order parameter, $\phi,$ is assumed to be dimensionless, as usual. With $x$ having units of length (we write $x\sim L$), the coefficients $c_{2j}$ have units of $L^{2j}$ as is clear from (2.3). The surface tension in the form (2.6) is defined as having units of length (as is typical in the reduced units used in physics) since

$c_{2j}\left\Vert D^{j}\phi\right\Vert _{L^{2}\left( \mathbb{R}\right) }% ^{2} = c_{2j}\int_{\mathbb{R}}\left\vert D^{j}\phi\right\vert ^{2}dx\ \sim L^{2j}L^{-2j+1} = L.$

The variables with units of length are $c_{2n}^{1/(2n)}$ and $\sigma$. We define

$z_{1}: = c_{2}^{1/2}, ..., z_{N}: = c_{2N}^{1/\left( 2N\right) }$

A set of dimensionless variables are

$\Pi_{0} = \frac{\sigma}{z_{1}}, \Pi_{1} = \frac{z_{2}}{z_{1}}% , ..., \Pi_{N-1} = \frac{z_{N}}{z_{1}}.$

(3.1)

The basic principle of dimensional analysis [1] is that a dimensionless quantity such as $\Pi_{0}$ can only depend on other dimensionless quantities, namely, $\Pi_{1}, ..., \Pi_{N-1}$ through some function $G:$

$\Pi_{0} = G\left( \Pi_{1}, ..., \Pi_{N-1}\right)$

(3.2)

Assuming a smooth solution $\phi$ to the phase field equation, (2.3), one can see from (2.5) or (2.6) that $G$ is differentiable in $z_{1}$ so we can write

$\frac{\partial\Pi_{0}}{\partial z_{1}} = % %TCIMACRO{\dsum \limits_{j = 1}^{N-1}}% %BeginExpansion {\displaystyle\sum\limits_{j = 1}^{N-1}} %EndExpansion \frac{\partial G}{\partial\Pi_{j}}\frac{\partial\Pi_{j}}{\partial z_{1}} = % %TCIMACRO{\dsum \limits_{j = 1}^{N-1}}% %BeginExpansion {\displaystyle\sum\limits_{j = 1}^{N-1}} %EndExpansion \frac{\partial G}{\partial\Pi_{j}}\left( -\frac{z_{j+1}}{z_{1}^{2}}\right) \ .$

(3.3)

For $k\geq2,$ we have similarly

$\frac{\partial\Pi_{0}}{\partial z_{k}} = \frac{\partial G}{\partial\Pi_{k-1}% }\frac{\partial\Pi_{k-1}}{\partial z_{k}} = \frac{\partial G}{\partial\Pi_{k-1}% }\left( \frac{1}{z_{1}}\right)$

(3.4)

Using (3.3) and (3.4) together, one obtains

$\frac{\partial\Pi_{0}}{\partial z_{1}} = % %TCIMACRO{\dsum \limits_{j = 1}^{N-1}}% %BeginExpansion {\displaystyle\sum\limits_{j = 1}^{N-1}} %EndExpansion z_{1}\frac{\partial\Pi_{0}}{\partial z_{j+1}}\left( -\frac{z_{j+1}}{z_{1}% ^{2}}\right) \ .$

(3.5)

At this point we regard $\sigma$ as a function of the $z_{j}$ rather than the $c_{2j}$. In other words, we define $\tilde \sigma \left ({{z_1}\%, ..., {z_N}} \right) = \sigma \left ({{c_2}, ..., {c_{2N}}} \right)$ and subsequently drop the tilda, as we will only use surface tension as a function of the $z_{j}$ below. Also, by using the definition of $\Pi_{0}$ we have for $k\geq2,$ we use the identity expressed by the definition of $\Pi_{0}$ so $\Pi_{0} =$ $\sigma\left (z_{1}, ..., z_{n}\right) /z_{1}$ with derivatives

$\begin{gathered} \frac{{\partial {\Pi _0}}}{{\partial {z_1}}} = \frac{{\partial \left( {\frac{\sigma }{{{z_1}}}} \right)}}{{\partial {z_1}}} = - \frac{1}{{z_1^2}}\sigma + \frac{1}{{{z_1}\% }}\frac{{\partial \sigma }}{{\partial {z_1}}}\;, \hfill \\ \frac{{\partial {\Pi _0}}}{{\partial {z_k}}} = \frac{{\partial \left( {\frac{\sigma }{{{z_1}}}} \right)}}{{\partial {z_k}}} = \frac{1}{{{z_1}}}\frac{{\partial \sigma }}{{\partial {z_k}}}\;. \hfill \\ \end{gathered}$

(3.6)

Thus one obtains

$\begin{gathered} - \frac{1}{{z_1^2}}\sigma + \frac{1}{{{z_1}}}\frac{{\partial \sigma }}{{\partial {z_1}}} = - \frac{1}{{{z_1}}}\sum\limits_{j = 1}^{N - 1} {{z_j}_{ + 1}} \frac{{\partial {\Pi _0}}}{{\partial {z_{j + 1}}}} = - \frac{1}{{{z_1}}}\sum\limits_{j = 1}^{N - 1} {{z_j}_{ + 1}} \frac{1}{{{z_1}}}\frac{{\partial \sigma }}{{\partial {z_{j + 1}}}} \hfill \\ = - \frac{1}{{z_1^2}}\sum\limits_{j = 2}^N {{z_j}} \frac{{\partial \sigma }}{{\partial {z_j}}}. \hfill \\ \end{gathered}$

Rewriting this by multiplying by $z_{1}^{2}$ we have

$%TCIMACRO{\dsum \limits_{j = 1}^{N}}% %BeginExpansion {\displaystyle\sum\limits_{j = 1}^{N}} %EndExpansion z_{j}\frac{\partial\sigma}{\partial z_{j}} = \sigma, \ \ or, \ z\cdot\nabla \sigma = \sigma$

(3.7)

Hence this is in the form of a Clairaut's equation [8].

---

4. Solution to Clairaut's equation for surface tension.

Clairaut's equation, (3.7) is a first order nonlinear partial differential equation that can be solved by the well-known method of characteristics (see e.g., [8]). In order to obtain a solution in $\mathbb{R}^{N}$ (actually the portion of $\mathbb{R}^{N}$ where each $z_{i}$ is non-negative) we will need to have a surface of "initial conditions" $\Gamma$ described below.

Definition. An acceptable manifold $\Gamma$ is a smooth surface in $\mathbb{R}^{N}$ such that the tangent plane to $\Gamma$ at any point $\vec {s}\in\Gamma$ intersects at a nonzero angle with any ray $\vec{v}$ emanating from the origin and going through an arbitrary point $\left ({z_1^0\%, z_2^0, ..., z_N^0} \right)$ such that $z_{i}^{0}\geq0$.

We define the coordinates $\left (r, s\right)$ with $r\in\mathbb{R}$, $s\in\mathbb{R}^{N-1}$ that form a new coordinate system

$\left( r\left( z_{1}, .., z_{N}\right), s\left( z_{1}, .., z_{N}\right) \right), \ \ \hat{\sigma}\left( r, s\right) = \sigma\left( z_{1}% , ..., z_{N}\right)$

(4.1)

and use the chain rule to obtain,

$\frac{\partial\hat{\sigma}}{\partial r} = \frac{\partial\sigma}{\partial z_{1}% }\frac{\partial z_{1}}{\partial r}+...+\frac{\partial\sigma}{\partial z_{N}% }\frac{\partial z_{N}}{\partial r}.$

One has then, from the standard methods of characteristics,

$\begin{gathered} \frac{{d{z_1}}}{{dr}} = {z_1}\;, \;so, \;{z_1}\left( {r;s} \right) = {C_1}\left( s \right){e^r}, \hfill \\ ... \hfill \\ \frac{{d{z_N}}}{{dr}} = {z_N}\;, \;so, \;{z_N}\left( {r;s} \right) = {C_N}\left( s \right){e^r}, \hfill \\ \frac{{d\sigma }}{{dr}} = \sigma, \;so, \;\sigma \left( {r;s} \right) = {C_{N + 1}}\left( s \right){e^r}. \hfill \\ \end{gathered}$

(4.3)

Next, we need to satisfy the "initial conditions" i.e., the values on the surface $\left (z_{1}, ..., z_{n}\right) \in\Gamma\subset\mathbb{R}^{N}$ where $\sigma$ is specified. For the $z_{i}$ variables, we make the choice that when $s\in\Gamma,$ we have $r = 0.$ Hence, we have $\sigma\left (0, s\right) = :\sigma_{0}\left (s\right)$ as the conditions on the surface $\Gamma,$ yielding the solutions

$\begin{gathered} {z_1}\left( {r;s} \right) = {C_1}\left( s \right){e^r}, ..., {z_N}\left( {r, s} \right) = {C_N}\left( s \right){e^r}, \hfill \\ \sigma \left( {r;s} \right) = {\sigma _0}\left( s \right){e^r}. \hfill \\ \end{gathered}$

with characteristics defined by

$\frac{z_{1}\left( r;s\right) }{z_{2}\left( r;s\right) } = \frac{C_{1}\left( s\right) }{C_{2}\left( s\right) }, ....$

This means that there is a fixed ratio of $z_{i}$ to $z_{j}$ for all $i, j\in\left\{ 1, ..., N\right\}$.

In summary, one can solve for all $\sigma\left (z_{1}, ..., z_{N}\right)$ in the subset $\mathbb{R}_{+}^{N}\subset\mathbb{R}^{N}$ for which all of the $z_{j}$ are nonnegative provided we specify the values of $\sigma\left (z_{1}, ..., z_{N}\right)$ on a surface $\Gamma\subset\mathbb{R}_{+}^{N}$ that is convex, and each point of $\Gamma$ intersects at a non-zero angle with each ray that emanates from the origin. If the surface $\Gamma$ is not convex, then the characteristics may intersect so that we obtain only local solutions.

---

5. Surface tension and homogeneous functions

We let $\Omega \subset\mathbb{R}_{+}^{N}$ be an open cone, i.e., if $z\in\Omega$ is in the set then so is $tz$, where $\mathbb{R}_ + ^N: = \left\{ {z \in \mathbb{R}{\% ^N}:{z_j} > 0} \right\}.$

Definition. For any scalar $k$ a real-valued function $f\left (z_{1}, ..., z_{n}\right)$ with $\left (z_{1}, ..., z_{n}\right) \in\Omega$ is homogeneous of degree $k$ if

$f\left( tz_{1}, ..., tz_{n}\right) = t^{k}f\left( z_{1}, ..., z_{n}\right) for\ all\ t>0.$

(5.1)

We recall two classical results.

Theorem. Let $f$ be a $C^{1}$ function on an open cone in $\mathbb{R}^{n}$. If $f$ is homogeneous of degree $k$ then its first order partial derivatives are homogeneous of degree $k-1.$

Theorem. Let $f:\Omega\rightarrow\mathbb{R}$ be a continuously differentiable function. Then the following are equivalent:

(1) $f$ is a homogeneous function of degee $k;$

(2) $f$ satisfies for $z\in\Omega,$ the equation

$\sum_{j = 1}^{n}z_{j}\frac{\partial f\left( z\right) }{\partial z_{j}% } = kf\left( z\right) .$

(5.2)

We now apply these concepts to the surface tension. Since $\sigma\left (z\right)$ satisfies (3.7), the theorem implies that $\sigma$ is homogeneous of degree 1. Hence, all of its partial derivatives are of degree 0. Hence we have the relation, for any $j\in\left\{ 1, ..., n\right\}$

$z_{1}\frac{\partial}{\partial z_{1}}\left( \frac{\partial\sigma}{\partial z_{j}}\right) +...+z_{n}\frac{\partial}{\partial z_{n}}\left( \frac {\partial\sigma}{\partial z_{j}}\right) = 0.$

(5.3)

Hence, each of the partial derivatives $\partial\sigma/\partial z_{j}$ satisfies the homogeneous Clairot's equation, $\vec{z}\cdot\vec{\nabla}% \sigma = 0$. Stated differently, the fact that the derivatives are homogeneous of order $0$ means that for any $t>0$ one has

$\frac{\partial\sigma}{\partial z_{j}}\left( tz_{1}, ..., tz_{n}\right) = \frac{\partial\sigma}{\partial z_{j}}\left( z_{1}, ..., z_{n}\right),$

i.e., the partial derivative is constant along the entire ray $\left\{ tz:t>0, z\in\mathbb{R}_{+}^{n}\right\}.$

One can also obtain similar results for a particular subspace. For example, if $n = 3,$ and we set $z_{2} = 0$ then one has similar results on the $z_{1}z_{3}$ plane.

---

Conflict of Interest

The author declares no conflicts of interest in this paper.

---