Interface dynamics for an Allen-Cahn-type equation governing a matrix-valued field

TL;DR

This paper analyzes the dynamics of a matrix-valued Allen-Cahn equation, revealing how solutions evolve towards orthogonal matrices and how interfaces move under curvature and surface diffusion effects.

Contribution

The paper introduces asymptotic methods to describe the solution behavior of a matrix-valued Allen-Cahn equation, including interface evolution driven by curvature and surface diffusion.

Findings

Solutions rapidly approach orthogonal matrices for single-signed determinant initial conditions.

Interface evolution is driven by curvature and surface diffusion, depending on matrix dimension.

Numerical experiments support the analytical asymptotic descriptions.

Abstract

We consider the initial value problem for the generalized Allen-Cahn equation, \[\partial_t \Phi = \Delta \Phi-\varepsilon^{-2} \Phi (\Phi^t \Phi - I), \qquad x \in \Omega, \ t\geq 0,\] where $\Phi$ is an $n\times n$ real matrix-valued field, $\Omega$ is a two-dimensional square with periodic boundary conditions, and $\varepsilon > 0$. This equation is the gradient flow for the energy, $E(\Phi) := \int \frac{1}{2} \|\nabla \Phi \|^2_F+\frac{1}{4 \varepsilon^2} \| \Phi^t \Phi - I \|^2_F$, where $\| \cdot \|_F$ denotes the Frobenius norm. The primary contribution of this paper is to use asymptotic methods to describe the solution of this initial value problem. If the initial condition has single-signed determinant, at each point of the domain, at a fast $O(\varepsilon^{-2} t)$ time scale, the solution evolves towards the closest orthogonal matrix. Then, at the $O(t)$ time scale, the solution evolves according to the $O_n$ diffusion equation. Stationary solutions to the $O_n$ diffusion equation are analyzed for $n=2$. If the initial condition has regions where the determinant is positive and negative, a free interface develops. Away from the interface, in each region, the matrix-valued field behaves as in the single-signed determinant case. At the $O(t)$ time scale, the interface evolves in the normal direction by curvature. At a slow $O(\varepsilon t)$ time scale, the interface is driven by curvature and the surface diffusion of the matrix-valued field. For $n=2$, the interface is driven by curvature and the jump in the squared tangental derivative of the phase across the interface. In particular, we emphasize that the interface when $n\geq 2$ is driven by surface diffusion, while for $n=1$, the original Allen--Cahn equation, the interface is only driven by mean curvature. A variety of numerical experiments are performed to verify, support, and illustrate our analytical results.

Equations

$${\partial }_t\Phi =\Delta \Phi -{\epsilon }^{-2}\Phi ({\Phi }^t\Phi -I),x\in \Omega ,t\ge 0,$$

$$\{\begin{array}{ll}& {\partial }_{t}A=\Delta A-{\epsilon }^{-2}A(A^{t}A-I),x\in \Omega ,t>0\\ & A(t=0,x)=A_0(x),\end{array}$$

$$E(A):={\int }_{\Omega }\frac{1}{2}\Vert \nabla A{\Vert }_F^2+{\epsilon }^{-2}W(A),\text{where }W(A):=\frac{1}{4}\Vert A^{t}A-I{\Vert }_F^2,$$

$$\{\begin{array}{ll}{A}_0(x)>0,& x\in {\Omega }_{+}\\ A_0(x)<0,& x\in {\Omega }_{-}\end{array},$$

$$O_n=S{O}_n\cup S{O}_n^{-}$$

$${\partial }_{t}B(t,x)=\frac{1}{2}\left(\Delta B(t,x)B^t(t,x)-B(t,x)\Delta B^t(t,x)\right)B(t,x)$$

$$B(0,x)={\Pi }_{O_n}{A}_0(x).$$

$$\{\begin{array}{ll}\det (A_0(x))>0,& x\in {\Omega }_{+}\\ \det (A_0(x))<0,& x\in {\Omega }_{-}\end{array},$$

$$A={\bar{A}}_0+\epsilon {\bar{A}}_1+{\epsilon }^2{\bar{A}}_2+o({\epsilon }^2)$$

$$\bar{A}({\bar{A}}^t\bar{A}-I)$$

$$+{\epsilon }^2({\bar{A}}_0{\bar{A}}_0^t{\bar{A}}_2+{\bar{A}}_0{\bar{A}}_2^t{\bar{A}}_0+{\bar{A}}_2{\bar{A}}_0^t{\bar{A}}_0+{\bar{A}}_0{\bar{A}}_1^t{\bar{A}}_1+{\bar{A}}_1{\bar{A}}_0^t{\bar{A}}_1+{\bar{A}}_1{\bar{A}}_1^t{\bar{A}}_0-{\bar{A}}_2)+o({\epsilon }^2).$$

$${\bar{A}}_i={\bar{A}}_i(x,t,{\tau }_1),i=1,2.$$

$${\partial }_{t}A={\epsilon }^{-2}{\partial }_{{\tau }_1}{\bar{A}}_0+{\epsilon }^{-1}{\partial }_{{\tau }_1}{\bar{A}}_1+{\partial }_{{\tau }_1}{\bar{A}}_2+{\partial }_t{\bar{A}}_0+o(1).$$

$${\epsilon }^{-2}{\partial }_{{\tau }_1}{\bar{A}}_0+{\epsilon }^{-1}{\partial }_{{\tau }_1}{\bar{A}}_1+o({\epsilon }^{-1})$$

$$={\epsilon }^{-2}{\bar{A}}_0({\bar{A}}_0^t{\bar{A}}_0-I)+{\epsilon }^{-1}({\bar{A}}_0{\bar{A}}_0^t{\bar{A}}_1+{\bar{A}}_0{\bar{A}}_1^t{\bar{A}}_0+{\bar{A}}_1{\bar{A}}_0^t{\bar{A}}_0-{\bar{A}}_1)+o({\epsilon }^{-1}).$$

$${\partial }_{{\tau }_1}{\bar{A}}_0={\bar{A}}_0({\bar{A}}_0^t{\bar{A}}_0-I).$$

$${\partial }_{t}B=B(B^{t}B-I)$$

$$B(t=0)=B_0$$

$$U({\partial }_t\Sigma )V^t=U(\Sigma {\Sigma }^t\Sigma -\Sigma )V^t.$$

$${\partial }_{{\tau }_1}{\bar{A}}_1={\bar{A}}_0{\bar{A}}_0^t{\bar{A}}_1+{\bar{A}}_0{\bar{A}}_1^t{\bar{A}}_0+{\bar{A}}_1{\bar{A}}_0^t{\bar{A}}_0-{\bar{A}}_1.$$

$${\bar{A}}_1({\tau }_1,x)=0$$

$${\partial }_t{\bar{A}}_0+o(1)$$

$$+({\bar{A}}_0{\bar{A}}_0^t{\bar{A}}_2+{\bar{A}}_0{\bar{A}}_2^t{\bar{A}}_0+{\bar{A}}_2{\bar{A}}_0^t{\bar{A}}_0+{\bar{A}}_0{\bar{A}}_1^t{\bar{A}}_1+{\bar{A}}_1{\bar{A}}_0^t{\bar{A}}_1+{\bar{A}}_1{\bar{A}}_1^t{\bar{A}}_0-{\bar{A}}_2+\Delta {\bar{A}}_0)+o(1).$$

$${\bar{A}}_0({\bar{A}}_0^t{\bar{A}}_0-I)=0.$$

$${\bar{A}}_0^t{\bar{A}}_0=I={\bar{A}}_0{\bar{A}}_0^t$$

$${\bar{A}}_0{\bar{A}}_0^t{\bar{A}}_1+{\bar{A}}_0{\bar{A}}_1^t{\bar{A}}_0+{\bar{A}}_1{\bar{A}}_0^t{\bar{A}}_0-{\bar{A}}_1=0.$$

$${\partial }_t{\bar{A}}_0$$

$$=\Delta {\bar{A}}_0-({\bar{A}}_2+{\bar{A}}_0{\bar{A}}_2^t{\bar{A}}_0),$$

$$0$$

$$={\bar{A}}_0^t(\Delta {\bar{A}}_0)-{\bar{A}}_0^t({\bar{A}}_2+{\bar{A}}_0{\bar{A}}_2^t{\bar{A}}_0)+(\Delta {\bar{A}}_0{)}^t{\bar{A}}_0-({\bar{A}}_2^t+{\bar{A}}_0^t{\bar{A}}_2{\bar{A}}_0^t){\bar{A}}_0.$$

$${\bar{A}}_0^t{\bar{A}}_2+{\bar{A}}_2^t{\bar{A}}_0=\frac{1}{2}\left({\bar{A}}_0^t(\Delta {\bar{A}}_0)+(\Delta {\bar{A}}_0{)}^t{\bar{A}}_0\right).$$

$${\partial }_t{\bar{A}}_0$$

$$=\frac{1}{2}\left(\Delta {\bar{A}}_0{\bar{A}}_0^t-{\bar{A}}_0(\Delta {\bar{A}}_0{)}^t\right){\bar{A}}_0$$

$${\partial }_{t}B(t,x)$$

Full text

[TABLE]

where $\Phi$ is an $n\times n$ real matrix-valued field, $\Omega$ is a two-dimensional square with periodic boundary conditions, and $\varepsilon>0$. This equation is the gradient flow for the energy, $E(\Phi):=\int\frac{1}{2}\|\nabla\Phi\|^{2}_{F}+\frac{1}{4\varepsilon^{2}}\|\Phi^{t}\Phi-I\|^{2}_{F}$, where $\|\cdot\|_{F}$ denotes the Frobenius norm. The primary contribution of this paper is to use asymptotic methods to describe the solution of this initial value problem.

If the initial condition has single-signed determinant, at each point of the domain, at a fast $O(\varepsilon^{-2}t)$ time scale, the solution evolves towards the closest orthogonal matrix. Then, at the $O(t)$ time scale, the solution evolves according to the $O_{n}$ diffusion equation. Stationary solutions to the $O_{n}$ diffusion equation are analyzed for $n=2$.

If the initial condition has regions where the determinant is positive and negative, a free interface develops. Away from the interface, in each region, the matrix-valued field behaves as in the single-signed determinant case. At the $O(t)$ time scale, the interface evolves in the normal direction by curvature. At a slow $O(\varepsilon t)$ time scale, the interface is driven by curvature and the surface diffusion of the matrix-valued field. For $n=2$, the interface is driven by curvature and the jump in the squared tangental derivative of the phase across the interface. In particular, we emphasize that the interface when $n\geq 2$ is driven by surface diffusion, while for $n=1$, the original Allen–Cahn equation, the interface is only driven by mean curvature.

A variety of numerical experiments are performed to verify, support, and illustrate our analytical results.

B. Osting is partially supported by NSF DMS 16-19755 and 17-52202.

X.-P. Wang was supported in part by the Hong Kong Research Grants Council (GRF grants 16302715, 16324416, 16303318, and NSFC-RGC joint research grant N-HKUST620/15)

1. introduction

We consider the initial value problem for the generalized Allen–Cahn equation,

[TABLE]

where $A(t,x)\in M(n)$ is a real matrix-valued field and $\varepsilon>0$ is a small parameter. For simplicity, we take the domain $\Omega$ to be a two-dimensional square, $[-1/2,1/2]^{2}$, with periodic boundary conditions. It is not difficult to show that (1) is the gradient flow for the energy,

[TABLE]

and $\|\cdot\|_{F}$ denotes the Frobenius norm. Roughly speaking, for $\varepsilon$ small, the solution to (1) is smoothed by the first term and the second term keeps the pointwise values of the matrix-valued field near $O_{n}$, the $n\times n$ orthogonal matrix group. This problem was first introduced in [OW17] as a model problem for several applications where a smooth matrix-valued field arises, including crystallography, where the matrix-valued field describes the local crystal orientation and inverse problems in image analysis e.g., diffusion tensor MRI or fiber tractography, where it is of interest to estimate a matrix- or orientation-valued function [OW18].

It is clear that when $n=1$, (1) reduces to the original Allen–Cahn equation [CH13], which models the behavior of two immiscible fluids. In (1), because the reaction rate is large compared to the diffusion rate, the solution at each point $x\in\Omega$ quickly tends to a stable equilibrium state of the reaction process, i.e., a minimum of $W(A)$. For the case $n=1$, the local minima of $W(A)$ are $1$ and $-1$. There are two cases.

• (i)

If $A_{0}(x)>0$ (or $A_{0}(x)<0$ respectively) for every $x\in\Omega$, then $A(t,x)$ will tend to $1$ (or $-1$ resp.). In this case, the effect of diffusion is only to slightly change the rate at which the solution approaches $1$ (or $-1$ resp.).

• (ii)

However, if $\Omega_{+}$ and $\Omega_{-}$ are such that $\Omega=\overline{\Omega_{+}}\cup\overline{\Omega_{-}}$ and $\Omega_{+}\cap\Omega_{-}=\emptyset$, with

[TABLE]

then a boundary layer in the solution develops at an interface between the two subdomains. Through a boundary layer expansion, one can show that $A(t,x)=\pm 1$ for $x$ away from the interface and that the interface evolves in the normal direction by its mean curvature.

We refer to [BK91, BR93] and references therein for more details on the $n=1$ case.

For $n\geq 2$, the minimizers of $W(A)$ are elements of $O_{n}$. Recall that

[TABLE]

where $SO_{n}$ denotes the special orthogonal group of the orthogonal matrices with determinant $1$ and $SO_{n}^{-}$ is the set of orthogonal matrices with determinant $-1$. In this paper, we use matched asymptotic expansion methods [E11] to show that, as in the $n=1$ case, there are two cases.

• (i)

If the initial condition $A_{0}(x)$ has either positive or negative determinant for each $x\in\Omega$, then no interface develops. We show in Section 2.2 that the $O(t)$ time dynamics of the leading order solution satisfies the $O_{n}$ diffusion equation,

[TABLE]

with initial condition given by

[TABLE]

Here, and throughout this paper, we use $\Pi_{O_{n}}A=\arg\min_{B\in O_{n}}\|A-B\|_{F}$ to denote the closest point in the orthogonal matrix group to the matrix $A$. We show in Proposition 2.2 that the leading order solution remains in $O_{n}$ pointwise for all time $t>0$. We discuss stationary solutions to the $O_{n}$ diffusion equation in Section 2.3.

• (ii)

In the second case, the initial condition $A_{0}(x)$ satisfies

[TABLE]

for $\Omega_{+}$ and $\Omega_{-}$ satisfying $\Omega=\overline{\Omega_{+}}\cup\overline{\Omega_{-}}$ and $\Omega_{+}\cap\Omega_{-}=\emptyset$. In this case, when $x$ is away from the interface, the behavior is similar to the first case. We derive a motion law for the interface at two time scales. At the $O(t)$ time scale, the interface evolves in the normal direction by curvature, as in the $n=1$ case; see Section 3.1. At a slow $O(\varepsilon t)$ time scale, we show in Section 3.2 that the interface is driven by the surface diffusion of the matrix-valued field and the curvature. For $n=2$, in Proposition 3.3, we show that the surface diffusion term can be written as the jump in the squared tangental derivative of the phase across the interface.

In particular, we emphasize that the interface when $n\geq 2$ is driven by surface diffusion, while for $n=1$, the original Allen–Cahn equation, the interface is only driven by mean curvature.

The results obtained via asymptotics are verified, supported, and illustrated in Section 4 through a wide variety of numerical experiments.

The model (1) considered in this paper can be viewed as a special case of the general model studied in [Lin+12]. In [Lin+12], an energy of the form in (2) is considered for high-dimensional, vector-valued functions and general assumptions on the potential $W$. General results for the phase transition of stationary solutions between minima of $W$ are derived. In the present paper, we consider time dynamics for our specific model.

Outline

The paper is organized as follows. In Section 2, we derive the behavior of the matrix-valued field satisfying (1) if the initial field, $A_{0}(x)$, only takes values in $SO_{n}$ or $SO_{n}^{-}$. In Section 3, we discuss the case when an interface develops between subdomains where $\det(A(x,t))>0$ and $\det(A(x,t))<0$. We develop a boundary layer around the interface and derive the motion of the interface at different time scales. Some numerical experiments are performed in Section 4. We conclude with a discussion in Section 5.

2. Evolution of an initial matrix-valued field with single-signed determinant

In this section, we discuss the case where the initial matrix-valued field, $A_{0}(x)$, is continuous and has positive determinant at each point $x\in\Omega$. The case where $A_{0}(x)$ has negative determinant everywhere is analogous. In this case, there is no interface appearing in the dynamics of the system. We consider the asymptotic expansion

[TABLE]

and the initial condition $A(t,x)|_{t=0}=A_{0}(x)$, which we assume to appear at the $O(1)$ scale. Then, we expand the nonlinear term on the right hand side of (1),

[TABLE]

We take two time scales: $t$ and a fast time scale $\tau_{1}:=\varepsilon^{-2}t$ and write

[TABLE]

We have $\partial_{t}=\varepsilon^{-2}\partial_{\tau_{1}}$ and

[TABLE]

We insert our ansatz into (1) and collect terms at each order in $\varepsilon$. Using (2) and (4) in (1) yields

[TABLE]

2.1. Behavior at the $O(\varepsilon^{-2}t)$ time scale

Collecting $O(\varepsilon^{-2})$ terms in (5) yields

[TABLE]

For $n=1$, it is well known that the solution of (6) approaches $1$ if the initial value is positive and approaches $-1$ if the initial value is negative as $\tau_{1}\rightarrow\infty$.

For $n\geq 2$, at each point $x\in\Omega$, as $\tau_{1}\rightarrow\infty$, the solution of (6) approaches a matrix in $SO_{n}$ if the initial matrix has a positive determinant and approaches a matrix in $SO_{n}^{-}$ if the initial matrix has a negative determinant, as shown in the following Lemma.

Lemma 2.1**.**

For the dynamic system

[TABLE]

for a non-singular initial $n\times n$ matrix $B_{0}$ ($n\geq 2$), as $t\rightarrow\infty$, the solution $B(t)$ approaches the nearest orthogonal matrix to $B_{0}$, written $\Pi_{O_{n}}B_{0}$.

Proof.

Write the singular-value decomposition of $B_{0}$ as $B_{0}=U\Sigma_{0}V^{t}$, where the diagonal values of $\Sigma_{0}$ are denoted by $\sigma_{0,i}$ ($i\in[n]$). Since the right hand side of the equation can initially be written as $U(\Sigma_{0}\Sigma_{0}^{t}\Sigma_{0}-\Sigma_{0})V^{t}$, the solution $B(t)$ also admits a singular-value decomposition with the same $U$ and $V$. Then, we can write the dynamic system as

[TABLE]

For each diagonal element $\sigma_{i}$ $(i\in[n])$ in $\Sigma$, we have $\partial_{t}\sigma_{i}=\sigma_{i}^{3}-\sigma_{i}$. Since the $\sigma_{i}(0)>0$ for $i\in[n]$, we have $\sigma_{i}(t)\rightarrow 1$ as $t\rightarrow\infty$. That implies, as $t\rightarrow\infty$, $B(t)$ approaches $UV^{t}$, which is the closest orthogonal matrix to $B_{0}$; see, e.g., [OW17, Lemma 1.1] ∎

Collecting $O(\varepsilon^{-1})$ terms in (5) yields

[TABLE]

Since $\bar{A}_{1}(0,x)=0$, we have

[TABLE]

is the solution to (7).

2.1.1. Summary of the behavior at the $O(\varepsilon^{-2}t)$ time scale

1. (1)

If the determinant of the initial matrix-valued field is positive for all $x\in\Omega$, the leading order matrix $\bar{A}_{0}(\tau_{1},x)$ at each point approaches the closest orthogonal matrix in $SO_{n}$. If the determinant of the initial matrix-valued field is negative for all $x\in\Omega$, the leading order matrix $\bar{A}_{0}(\tau_{1},x)$ at each point will approach the closest orthogonal matrix in $SO_{n}^{-}$. 2. (2)

The second order matrix $\bar{A}_{1}(\tau_{1},x)$ is [math] for any $\tau_{1}\geq 0$ and $x\in\Omega$.

2.2. Behavior at the $O(t)$ time scale

Using (2) and (4) in (1), we have at the time scale $O(t)$,

[TABLE]

Collecting $O(\varepsilon^{-2})$ terms in (8) yields

[TABLE]

Since $\bar{A}_{0}$ is non-singular, this implies that

[TABLE]

This means that at each point $x\in\Omega$, for any non-singular initial matrix field, the leading order $\bar{A}_{0}$ immediately approaches a matrix field with values in $O_{n}$. This can be interpreted as the long time behavior of the dynamics at the time scale $O(\varepsilon^{-2}t)$.

Collecting $O(\varepsilon^{-1})$ terms in (8), we obtain

[TABLE]

This is also consistent with the solution $\bar{A}_{1}(t,x)=0$ at the $O(\varepsilon^{-2}t)$ time scale.

At $O(1)$ in (8), we obtain

[TABLE]

where we have used the fact $\bar{A}_{1}=0$. We now take the derivative of (9) and use (10b) to obtain

[TABLE]

Using (9) and rearranging, we obtain

[TABLE]

Multiplying on the left by $\bar{A}_{0}$ and using (9), we insert this back into (10b) to obtain

[TABLE]

The following proposition shows that if initially $\partial_{t}\bar{A}_{0}$ is in $SO_{n}$ pointwise, then it will remain there for all time $t>0$.

Proposition 2.2**.**

We consider the initial value problem for the $O_{n}$ diffusion equation,

[TABLE]

where $B_{0}\colon\Omega\to O_{n}$ is given. Then $B(t,x)\in O_{n}$ for all $t\geq 0$.

Proof.

We compute

[TABLE]

∎

Remark 2.3*.*

We refer to (11) and (12) as the $O_{n}$ diffusion equation because it can be obtained from the diffusion equation for a matrix-valued field when the matrix is constrained to be $O_{n}$-valued. That is, introducing a Lagrange multiplier for the constraint and then solving for it yields precisely this equation.

Remark 2.4*.*

For $n=2$ and the ansatz,

[TABLE]

we compute

[TABLE]

We conclude that $B(x,t)$ satisfies the orthogonal diffusion equation (12) if $\eta_{t}=\Delta\eta$. The spherical diffusion equation is given by

[TABLE]

see, e.g., [EW00]. Making the ansatz $\phi(x,t)=e^{i\eta(t,x)}$, we find that $\eta_{t}=\Delta\eta$. Thus, we conclude that the $n=2$ orthogonal diffusion equation (12) with initial condition taking values in $SO_{n}$ is equivalent to the spherical diffusion equation. Due to this connection, we refer to $\eta=\eta(x,t)$ as the phase of the matrix-valued field, $B(x,t)$.

2.2.1. Summary of the behavior at the $O(t)$ time scale:

1. (1)

The leading order solution $\bar{A}_{0}$ take values in the orthogonal matrix group for all time. 2. (2)

The second order solution is $\bar{A}_{1}=0$. 3. (3)

The time dynamics of $\bar{A}_{0}$ is governed by the $O_{n}$ diffusion equation,

[TABLE]

Here the initial condition is pointwise the closest point to $A_{0}(x)$ in $O_{n}$.

2.3. Harmonic orthogonal matrix-valued fields

In this section, we consider stationary solutions of the $O_{n}$ diffusion equation (12), which we refer to as harmonic orthogonal matrix-valued fields, satisfying

[TABLE]

Note that (13a) just states that $(\Delta B)^{t}B$ is a symmetric matrix.

For $n=1$, the only solutions are $B=\pm 1$.

For $n=2$, for unknown phase $\eta\colon\Omega\to\mathbb{R}$, we consider the $SO_{n}$ ansatz,

[TABLE]

We compute

[TABLE]

We observe that $(\Delta B)^{t}B$ is symmetric if and only if $\Delta\eta=0$. We conclude that there exists a family of harmonic $SO_{n}$-valued fields on the torus of the form (14) where the phase is given by

[TABLE]

Several numerical experiments are performed in Section 4.2.1 to show that such fields are stationary for (1) and to investigate what happens if perturbations of such fields are taken as initial conditions.

3. Evolution of an $O_{n}$-valued initial field

We consider an initial condition $A_{0}(x)$ that satisfies

[TABLE]

for $\Omega_{+}$ and $\Omega_{-}$ satisfying $\Omega=\overline{\Omega_{+}}\cup\overline{\Omega_{-}}$ and $\Omega_{+}\cap\Omega_{-}=\emptyset$.

Denote $\Gamma=\bar{\Omega}_{+}\cap\bar{\Omega}_{-}$ and assume $\Gamma(t)$ is a finite collection of simple, closed, smooth curves in $\mathbb{R}^{2}$ so that we can find a parametric representation, at least locally, of the form

[TABLE]

We assume $s$ to be the arc-length parameter so that we have

[TABLE]

where $\vec{T}$ denotes the unit tangent vector, $\vec{n}$ denotes the unit outer normal vector, and $\kappa$ denotes the curvature (see Figure 1 for a diagram on the unit outer normal vector).

We introduce local coordinates near $\Gamma$ as follows. We assume that for every point $x$ in a neighborhood of $\Gamma$, there is a unique point $\vec{\varphi}(s,t)$ which is the orthogonal projection of $x$ onto $\Gamma(t)$. We then define a unique normal signed distance $\rho(x,t)$ from $x$ to $\Gamma(t)$, $\rho(x,t)=(x-\vec{\varphi})\cdot\vec{n}$. We have a transformation from $(x,t)$ to $(s,r,t)$ defined by

[TABLE]

where $r=\rho(x,t)$. We summarize several identities for the transformation from $(x,t)$ to $(s,r,t)$ in the following Lemma. A proof of this Lemma can be found in [DD12] .

Lemma 3.1**.**

*For the transformation rule defined in (17), we have the following equalities:

1.The normal velocity of $\Gamma$ at $\vec{\varphi}(s,t)$ is given by $\vec{v}\cdot\vec{n}=-\frac{\partial\rho}{\partial t}$.

2. $\nabla_{x}s=\frac{1}{1+r\kappa}\vec{T}$, $\Delta_{x}s=-r\frac{\partial\kappa}{\partial s}\frac{1}{(1+r\kappa)^{3}}$.

3. $\nabla_{x}r=\vec{n}$, $\Delta_{x}r=\frac{\kappa}{1+r\kappa}$.

4. For any function $u(x,t)=\tilde{u}(s,r,t)$,*

[TABLE]

Below, we study the inner layer expansion to study the behavior around the interface at the time scales $O(t)$ and $O(\varepsilon t)$.

3.1. Behavior at the $O(t)$ time scale.

At the $O(t)$ time scale, when $x$ is away from the interface $\Gamma$, the behavior exactly reduces to the case studied in Section 2.2. That is corresponding to the outer layer expansion for the system (1). We don’t repeat the calculation and refer the results to the summary in Section 2.2.1.

The inner expansion requires rescaling the normal coordinate by $z=\frac{r}{\varepsilon}$. Assume the expansion of $\vec{\varphi}(s,t)$ and $\kappa$ are

[TABLE]

Writing $\Delta_{x}$ in terms of $(s,z)$ using (19) yields

[TABLE]

Consider the expansion

[TABLE]

and write

[TABLE]

Similar to (2), we have

[TABLE]

Substituting (20), (21), (22), and (3.1) into (1) yields

[TABLE]

Collecting the $O(\varepsilon^{-2})$ terms in (24), we obtain

[TABLE]

Matching the outer expansion gives the boundary conditions

[TABLE]

where $\bar{A}_{0}\in O_{n}$. Note that (25) is independent of $s$ and thus we solve (25) for each $s$ independently. For $n=2$, the following proposition gives the explicit solution to (25) with the boundary conditions in (26).

Proposition 3.2**.**

The solution to the second-order differential equation for the $2\times 2$ matrix field, $B\colon\mathbb{R}\to M(2)$

[TABLE]

is given by

[TABLE]

where $\xi_{1}=\frac{\eta_{-}+\eta_{+}}{2}$ and $\xi_{2}=\frac{\eta_{-}-\eta_{+}}{2}$.

Proof.

The proposed solution takes the form $B(z)=U_{1}DU_{2}^{t}$, where $U_{i}$ for $i=1,2$ are matrices in $SO_{2}$ that are independent of $z$ and $D=\begin{bmatrix}1&0\\ 0&\tanh\left(\frac{z}{\sqrt{2}}\right)\end{bmatrix}$ is a diagonal matrix. Since $\sigma(z)=\tanh\left(\frac{z}{\sqrt{2}}\right)$ satisfies $\partial_{zz}\sigma=\sigma^{3}-\sigma$, we have that $B(z)$ satisfies the differential equation. Using the two trigonometric identities,

[TABLE]

$B(z)$ satisfies the boundary conditions as $z\to\pm\infty$. ∎

Hence, for $n=2$, by Proposition 3.2, we explicitly get the solution for (25) coupled with boundary conditions (26),

[TABLE]

where $\xi_{1}=\xi_{1}(s,t)$ and $\xi_{2}=\xi_{2}(s,t)$ are determined from the phase of the outer solution, $\bar{A}_{0}(x,t)$, for each $s$ and $t$.

Collecting the $O(\varepsilon^{-1})$ terms in (24) yields

[TABLE]

Taking the Frobenius inner product with $\partial_{z}\tilde{A}_{0}$ on both sides of (28) yields

[TABLE]

We integrate the above equation with respect to $z$ from $-\infty$ to $\infty$. Integrating by parts, we can rewrite the first term on the right hand side as

[TABLE]

where the boundary terms vanish because $\partial_{z}\tilde{A}_{0}=0$ at $z=\pm\infty$ from the solution to the leading order expansion and $\bar{A}_{1}=0$ in the outer layer which corresponds to $\tilde{A}_{1}=0$ at $z=\pm\infty$ from the asymptotic matching. The second term on the right hand side can be rewritten as

[TABLE]

where the boundary terms vanish because of $\bar{A}_{1}=0$ in the outer layer. We then obtain

[TABLE]

Denote

[TABLE]

We note that when $n=1$, $\gamma$ is the surface tension on the interface between two different phases [Peg89]. Since we assume $\tilde{A}_{0}\in SO_{n}$ on one side of $\Gamma(t)$ and $\tilde{A}_{0}\in SO_{n}^{-}$ on the other, we have $\gamma(s)>0$. Using (25), the right hand side of (29) vanishes and we have the normal velocity is given by

[TABLE]

It follows that the interface evolves according to the mean curvature flow along its normal direction at the $O(t)$ time scale.

3.1.1. Summary of the behavior near the interface at the $O(t)$ time scale:

1. (1)

The leading order solution, $\tilde{A}_{0}$, transitions from a matrix in $SO_{n}^{-}$ to a matrix in $SO_{n}$ in the boundary layer of the interface. It satisfies (25) with the boundary conditions in (26). For $n=2$, the leading order solution is explicitly given by (27). 2. (2)

The interface moves in the normal direction by the leading order of curvature, i.e., $\vec{v}\cdot\vec{n}=-\kappa_{0}$.

3.2. Behavior at the $O(\varepsilon t)$ time scale.

Now, we study the behavior at the time scale $O(\tau_{2})$ where $\tau_{2}=\varepsilon t$. Then, we have $\partial_{t}=\varepsilon\partial_{\tau_{2}}$.

First, we consider the behavior when $x$ is away from the interface $\Gamma$. Consider the expansion

[TABLE]

and insert it into (1) to obtain

[TABLE]

Collecting the terms at different orders of $\varepsilon$ and using the behavior at the $O(t)$ time scale, we obtain that at the $O(\tau_{2})$ time scale,

[TABLE]

For $n=2$, as shown in Section 2.3, there is a $O_{2}$ harmonic leading order matrix field, $\bar{A}_{0}$ of the form

[TABLE]

where the phase $\eta(x)$ satisfies $\Delta\eta=0$. Here, the signs in the second column are chosen depending on whether $x\in\Omega_{+}$ or $x\in\Omega_{-}$.

To study the behavior near the interface, we consider the expansion

[TABLE]

and rewrite (22) using $\tau_{2}=\varepsilon t$ to obtain

[TABLE]

Substituting (20), (3.1), (35), and (34) into (1), we have

[TABLE]

Collecting the $O(\varepsilon^{-2})$ terms in (36) yields

[TABLE]

which is same as (25). The boundary conditions at $z=\pm\infty$ are as in (26). Hence, $\tilde{A}_{0}(s,z,\tau_{2})$ has the same transition profile as obtained from (27) in the boundary layer.

At $O(\varepsilon^{-1})$ in (36), we have

[TABLE]

Taking the Frobenius inner product with $\partial_{z}\tilde{A}_{0}$ on both sides of (28) and integrating both sides with respect to $z$ from $-\infty$ to $\infty$ yields:

[TABLE]

Using (37) again leads us to

[TABLE]

Inserting (39) back to (38) yields

[TABLE]

Coupling with the boundary conditions from the outer expansion: $\tilde{A}_{1}(s,z,\tau_{2})=0$ at $z=\pm\infty$ implies that

[TABLE]

Collecting the $O(1)$ terms in (36), we obtain

[TABLE]

Using (39) and (40) simplifies (3.2) to

[TABLE]

Taking the Frobenius inner product with $\partial_{z}\tilde{A}_{0}$ on both sides of (42) and integrating both sides with respect to $z$ from $-\infty$ to $\infty$ yields:

[TABLE]

Combining with (37) and using the definition of $\gamma$ in (30) leads us to

[TABLE]

which gives the motion law for the interface at the $O(\tau_{2})$ time scale.

Proposition 3.3**.**

For $n=2$, $\tilde{A}_{0}$ admits the transition profile in (27) where $\xi_{1}(s,t)=\frac{\eta_{-}(s,t)+\eta_{+}(s,t)}{2}$ and $\xi_{2}(s,t)=\frac{\eta_{-}(s,t)-\eta_{+}(s,t)}{2}$. Here, $\eta_{+}(s,t)$ and $\eta_{-}(s,t)$ are determined from the outer solution, $\bar{A}_{0}(x,t)$, as in (33). The motion law in (43) simplifies to

[TABLE]

where $\left[\eta_{s}^{2}\right]_{\Gamma}=\left(\left(\partial_{s}\eta_{+}(s,t)\right)^{2}-\left(\partial_{s}\eta_{-}(s,t)\right)^{2}\right)$ is the jump in the squared tangental derivative of the phase across the interface $\Gamma$ and $\bar{\gamma}=\int_{-\infty}^{\infty}\left(1-\tanh^{2}(z/\sqrt{2})\right)^{2}\ dz$.

Proof.

In the transition profile (27), we write

[TABLE]

Here we write $\sigma(z)=\tanh\left(\frac{z}{\sqrt{2}}\right)$. We compute

[TABLE]

and

[TABLE]

Then we compute

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Using

[TABLE]

along with the above identities, we have

[TABLE]

and

[TABLE]

which is independent of $s$ and denoted by $\bar{\gamma}$. ∎

3.2.1. Summary of the behavior at the $O(\varepsilon t)$ time scale:

1. (1)

Away from the interface, $\bar{A}_{0}$, $\bar{A}_{1}$, and $\bar{A}_{2}$ satisfy (32). In particular, when $n=2$, $\bar{A}_{0}$ takes the form in (33) where the phase $\eta$ satisfies $\Delta\eta=0$. 2. (2)

The interface moves in the normal direction according to the motion law given in (43). In the case where $n=2$, the second term of the motion law is the jump in the squared tangental derivative of the phase across the interface $\Gamma$ as in (44).

4. Numerical experiments

In this section, we perform a variety of numerical experiments to support, verify, and illustrate our analytical results in Section 2 and Section 3. The algorithm we use is summarized in Section 4.1 and the numerical examples are described in Section 4.2.

4.1. Algorithm to solve (1) and implementation details.

To numerically solve (1), we use an efficient diffusion generated method recently developed in [OW17]. This method generalizes the Merriman-Bence-Osher method for mean curvature flow [Mer+94] and methods for the Ginzburg-Landau energy [Ruu+01, VO19]. The algorithm alternates a diffusion and a projection step as summarized in Algorithm 1. In [OW17], the Lyapunov function of Esedoglu and Otto [EO15] was extended to show that the method is non-increasing on iterates and hence, unconditionally stable. It was also proven that the spatially discretized iterates converge to a stationary solution in a finite number of iterations. We refer to [OW17] for more details and properties of the algorithm.

We implemented the algorithm in MATLAB. In all experiments, we consider the case when $n=2$ on a flat torus $\Omega=[-1/2,1/2]^{2}$ discretized using $1024\times 1024$ uniform grid points and set $\tau=0.015625$. The heat diffusion equation in Algorithm 1 is efficiently solved using the fast Fourier transform (FFT). The convergence criteria of the algorithm is taken to be

[TABLE]

for $tol=10^{-6}$. All reported results were obtained on a laptop with a 2.7GHz Intel Core i5 processor and 8GB of RAM.

Here we visualize an $O_{2}$ valued field by plotting the vector field generated by the first column vector. The second column vector is orthogonal to the first and the direction is indicated by color, when necessary.

4.2. Numerical examples

4.2.1. Evolution of $SO_{n}$-valued fields

We first perform a numerical experiment to verify the results in Section 2 for the time evolution of a single-signed determinant initial matrix-valued field. Without loss of generality, we consider the case where the initial matrix-valued field takes values in $SO_{n}$.

Figures 2 and 3 display the evolution of an $SO_{2}$ matrix-valued field with the initial condition given by

[TABLE]

for different choices of $\eta\colon\Omega\to\mathbb{R}$.

In Figure 2, we take

[TABLE]

From Figure 2, we see that the matrix-valued field evolves toward a uniform matrix-valued field, which, as discussed in Section 2.3, is a stationary state of the $O_{n}$ diffusion equation.

In Figure 3, we set

[TABLE]

We observe that the field evolves toward a field with $\eta(x)=2\pi x_{1}$. Again, since $\Delta\eta=0$, this is a harmonic orthogonal matrix-valued field; see Section 2.3.

To better understand the behavior in Figures 2 and 3, we recall the definition of the pair of indices of a matrix-valued field discussed in [OW17]. Let $v\colon\Omega\to\mathbb{C}$ be a complex-valued field with no zeros. Let $\gamma\colon[0,1]\rightarrow\Omega$ be a closed curve. We define the index of $\gamma$ with respect to $v$ to be

[TABLE]

Clearly the index of $\gamma$ is an integer and varies continuously with deformations to $\gamma$, so it depends only on the homotopy class of $\gamma$. For a torus, we can parameterize the homotopy classes by the number of times the curve wraps around $\Omega$ in the $x_{1}$- and $x_{2}$-directions. Furthermore, if we let $[\gamma]_{m,n}$ denote the equivalence class of curves that wraps around $\Omega$ $m$ times in the $x_{1}$-direction and $n$ times in the $x_{2}$-direction, then it is not difficult to see that

[TABLE]

So we can characterize the index of any curve in terms of the indices of $[\gamma]_{1,0}$ and $[\gamma]_{0,1}$. For a given field $v$, we let

[TABLE]

be the index pair corresponding to curves that wrap around $\Omega$ once in the $x_{1}$- and $x_{2}$-directions. For a matrix-valued field $A\colon\Omega\to SO_{2}$ or $A\colon\Omega\to SO_{2}^{-}$, we define the index pair, $I$, to be the index pair for the first column of $A$. For example, for the harmonic orthogonal matrix fields in (15), the index pair is $I=(n_{1},n_{2})$.

In Figure 2, the index pair for the initial condition is $(0,0)$ and the field evolves toward the harmonic orthogonal matrix field with index pair $(0,0)$, the uniform matrix field. In Figure 3, the index pair for the initial condition is $(1,0)$ and the field evolves toward the harmonic orthogonal matrix field with index pair $(1,0)$. We observe and generally expect that the index pair is invariant under flow by the $O_{n}$ diffusion equation.

4.2.2. Evolution of $O_{n}$-valued fields at the $O(t)$ time scale.

In this section, we check the motion law we derived in Section 3.1. That is, at the $O(t)$ time scale, if there is a line defect initially, the motion of the interface is driven by the curvature at each point. Note that at this time scale, we don’t see the effect from the matrix-valued field on the motion law of the interface. So we perform two experiments where the initial condition has the same line defect, but different initial matrix-valued fields. Specifically, we choose the following initial condition for different choices of $\eta\colon\Omega\to\mathbb{R}$,

[TABLE]

where $(r,\theta)$ is the corresponding polar coordinate of $x=(x_{1},x_{2})\in\Omega$.

In all subsequent figures, the domain is colored by the sign of the determinant of the matrix. For a matrix field $A\in H^{1}\left(\Omega;O_{n}\right)$, we use the convention

[TABLE]

In Figure 4, we display several snapshots of the time evolution for two different initial conditions. In the first column of Figure 4, the initial field is chosen as in (46) with $\eta(x)=\frac{\pi}{2}\sin(2\pi x_{1})$ and, in the second column, the initial field has $\eta(x)=2\pi x_{1}$. Hence the pair of indices of the initial field in the first column is $(0,0)$ and the pair of indices of the initial field in the second column is $(1,0)$. In both columns of Figure 4, we observe that the region where $A\in SO_{n}$ shrinks with the interface becoming a circle before vanishing. We observe that the time dynamics of the line defect for the two different initial conditions are very close. This is consistent with our analytical results in Section 3.1, that is, at the $O(t)$ time scale, the motion law is the leading order of the curvature of the line defect, which is independent of the matrix-valued field. For the evolution in the left column of Figure 4, the field continues to evolve toward the uniform solution for longer times than shown in the Figure.

4.2.3. Evolution of $O_{n}$-valued fields at the $O(\varepsilon t)$ time scale.

In this section, we check the motion law we derived in Section 3.2 at the $O(\varepsilon t)$ time scale. Note that, at the $O(\varepsilon t)$ time scale, we have the leading order of the curvature of the line defect satisfies $\kappa_{0}=0$ and we have $\Delta\eta=0$. Hence, we perform several experiments where the initial condition has two straight parallel line defects. Specifically, we choose the following initial condition

[TABLE]

for different choices of $\eta_{1},\eta_{2}\colon\Omega\rightarrow\mathbb{R}$ satisfying $\Delta\eta_{1}=0$ and $\Delta\eta_{2}=0$.

We first choose $\eta_{1}$ and $\eta_{2}$ so that $\eta_{1s}^{2}=\eta_{2s}^{2}$; the parallel line defects are stationary according to our analytical results in (44) in Section 3.2. In Figure 5, from the left to the right, we set

[TABLE]

Indeed, in Figure 5, all four columns show that the parallel line defects are stationary.

Figures 6 and 7 display snapshots of the time dynamics of the interface at different times for two choices of the phases $\eta_{1}$ and $\eta_{2}$. In Figure 6, we set the initial phases to be $\eta_{1}=2\pi x_{1}$ and $\eta_{2}=4\pi x_{1}$. In Figure 7, we set the initial phases to be $\eta_{1}=2\pi x_{1}$ and $\eta_{2}=8\pi x_{1}$. Thus, we have $[\eta_{s}^{2}]_{\Gamma}=-12\pi^{2}$ in Figure 6 and $[\eta_{s}^{2}]_{\Gamma}=-60\pi^{2}$ in Figure 7. In both figures, we observe that the straight line defects have nonzero speed along their normal directions. Comparing Figure 6 and Figure 7, we observe that the speed of the line defect in Figure 7 is about five times the speed of the line defect in Figure 6. All these observations are consistent with our analytical results on the motion law (44) in Section 4.2.3.

Figure 8 displays the snapshots of the time dynamics of the interface at different times where the initial phases are give by $\eta_{1}=8\pi x_{1}$ and $\eta_{2}=2\pi x_{1}$. We observe the dynamics in Figure 8 has opposite direction than in Figure 7. This is also consistent with our analytical result on the motion law (44) in Section 4.2.3.

5. Conclusion and discussion

In this paper, we used asymptotic methods to study the initial value problem for the generalized Allen-Cahn equation in (1). If the initial condition has single-signed determinant, at each point of the domain, at a fast $O(\varepsilon^{-2}t)$ time scale, the solution evolves towards the closest orthogonal matrix. Then, at the $O(t)$ time scale, the solution evolves according to the $O_{n}$ diffusion equation (12a). Stationary solutions to the $O_{n}$ diffusion equation were analyzed for $n=2$ in Section 2.3. If the initial condition has regions where the determinant is positive and negative, an interface develops. Away from the interface, in each region, the matrix-valued field behaves as in the single-signed determinant case. At the $O(t)$ time scale, the interface evolves in the normal direction by curvature. At a slow $O(\varepsilon t)$ time scale, for $n=2$, the interface is driven by curvature and the jump in the squared tangental derivative of the phase across the interface. In Section 4, we conducted a variety of numerical experiments to verify, support, and illustrate our analytical results.

In this paper, we have focused on the two-dimensional problem. We expect that the asymptotic methods in [DD12, DD14] could be used to study higher-dimensional problems. In this paper, we also only focused on a square with periodic boundary conditions. We used this to derive harmonic orthogonal matrix-valued-fields in Section 2.3 and in the numerical examples in Section 4. However, the asymptotic results from Sections 2 and 3 apply to other boundary conditions as well.

There are several places where we focused on the $n=2$ case. In particular, in Section 2.3, we derived explicit $O_{2}$ harmonic fields; in Proposition 3.2 we explicitly derived the transition profile for the boundary layer; and in Proposition 3.3, we were able to simplify the expression for the motion law at the slow $O(\varepsilon t)$ time scale in terms of the jump in the squared tangental derivative of the phase across the interface. It would be interesting to extend these results to $n\geq 3$.

Here, we have considered the $L^{2}$ gradient flow (1) of the energy $E$ in (2). An interesting equation would arise from considering the $H^{-1}$ gradient flow of $E$, a generalization of the Cahn–Hilliard equation [Peg89, DD12, DD14, CH13, Che+14, Wan+17]. Another way to generalize (1) would be to consider multi-phase systems as in [Rub+89, BR93].