# SPEECH-BASED EMOTION RECOGNITION WITH SELF-SUPERVISED MODELS USING ATTENTIVE
CHANNEL-WISE CORRELATIONS AND LABEL SMOOTHING

###### Abstract

When recognizing emotions from speech, we encounter two common problems: how
to optimally capture emotion-relevant information from the speech signal and
how to best quantify or categorize the noisy subjective emotion labels. Self-
supervised pre-trained representations can robustly capture information from
speech enabling state-of-the-art results in many downstream tasks including
emotion recognition. However, better ways of aggregating the information
across time need to be considered as the relevant emotion information is
likely to appear piecewise and not uniformly across the signal. For the
labels, we need to take into account that there is a substantial degree of
noise that comes from the subjective human annotations. In this paper, we
propose a novel approach to attentive pooling based on correlations between
the representations’ coefficients combined with label smoothing, a method
aiming to reduce the confidence of the classifier on the training labels. We
evaluate our proposed approach on the benchmark dataset IEMOCAP, and
demonstrate high performance surpassing that in the literature. The code to
reproduce the results is available at
github.com/skakouros/s3prl_attentive_correlation.

Index Terms—  emotion recognition, self-supervised features, iemocap, hubert,
wavlm, wav2vec 2.0

##  1 INTRODUCTION

Emotional expressions are a fundamental component of spoken interaction. When
we communicate with other people, we are implicitly monitoring their emotional
state and respond based on that emotional state [1]. Emotions fall in the
realm of prosodic function. In recent years, the importance of prosodic
qualities in speech has attracted increasing attention. This is also the case
for speech emotion recognition (SER) which has seen a growing interest with
the increasing role of spoken language interfaces in human-computer
interaction (HCI) applications [2]. However, recognizing emotions in speech
remains a challenging problem complicated by numerous factors including
fundamental issues of how emotion is defined, elicited, expressed, and
communicated [3] and extending to how we can capture this information from
speech.

The challenges in SER can be split into three distinct problems. First, we
encounter the issue of developing engineered representations that can robustly
capture the acoustic information in speech that best describe the variation
found across different emotions. This has been traditionally done using
features such as mel-frequency cepstral coefficients (MFCCs), filterbanks,
fundamental frequency, energy, zero-crossing rate, chroma-based features and
their feature functionals [4] or through standardized feature collections and
their functionals such as eGeMAPS [5]. More recently, self-supervised learning
(SSL) has shown its effectiveness in various domains, including SER, and is
becoming the new principle for extracting representations from speech. HuBERT
[6], Wav2vec 2.0 [7], WavLM [8], are some of the self-supervised approaches
for speech representation learning that have been used in the context of SER
[9, 5].

The second issue that we face in SER is the effective modeling of the long
temporal context over which emotions take place. Emotion specific information
lies beyond segmental productions and in longer time scales. These may include
parts of an utterance but can also span across one or more utterances. To
appropriate model long-term dependencies by capturing and connecting the
relevant cues across time suitable methods are necessary. These vary from
approaches that simply take the first and second order statistics of self-
supervised representations across time [9] to approaches that focus on complex
sequence modeling tasks [10]. For example Sarma et al. [10] used a TDNN
architecture combined with LSTM and self-attention to model the long-term
temporal context and to capture the emotionally relevant portions of speech.
In a recent work, Liu et al. [11] used a cascaded attention network to locate
the relevant emotion regions from the input features. Other approaches also
use different types of recurrent neural networks (RNNs) to explain the long
temporal contexts of emotions in speech [12].

The third and final issue comes from the observation that human emotional
expressions are often unclear and ambiguous, leading to disagreement and
confusion among human evaluators [13, 14]. This confusion might be partly
attributed to the multimodal nature of emotion expression. Facial expressions,
hand gestures, and speech with its prosodic and linguistic content all work
together in eliciting different emotions. Perhaps the absence of multimodality
may be one source of confusion that leads to overlaps in the clusters of the
different emotion classes. However, speech alone holds much relevant
information in its prosodic content that can be used for robust SER. Different
ways to tackle the problem of noisy labels and consequently the uncertainty in
predicting emotions have been suggested in the literature. These typically
include custom loss functions [15, 11] and modifications to the target hard
labels [16, 17]. For example Liu and Wang used a triplet loss to make anchor
utterances more similar to all other positive utterances [15] while Tarantino
et al. used regression targets instead of hard categorical targets by taking
the proportion of the classes within the annotations [17].

This paper presents a framework for SER that uses pre-trained speech models
with a novel approach to attentive pooling based on channel-wise correlations
on soft targets. We evaluate the framework with HuBERT [6], Wav2vec 2.0 [7],
and WavLM [8] upstream models. We use the SUPERB [9] evaluation setup
throughout our experiments. The effectiveness of our proposed framework is
evaluated on the interactive emotional dyadic motion capture (IEMOCAP) dataset
[18] and shows state-of-the-art performance in SER.

##  2 RELATION TO PRIOR WORK

The idea of taking the correlations between different filter responses over
the spatial extent of the feature maps to obtain a representation of the style
of an input image was introduced by Gatys et al. [19]. Their method was later
adapted to speech where it has found applications in speech generation and
voice conversion [20], pooling to obtain speaker embeddings [21], and
sentence-level tasks such as speaker identification, speaker verification, and
SER [22].

In this work, we extend the method for correlation pooling presented in [22].
In [22], it was shown that channel-wise correlations provide an alternative
way of extracting speaker and emotion information from self-supervised models,
providing also improvements over the standard mean and mean-std pooling (std
stands for standard deviation). In the present work, we add an attention
mechanism to pool representations before estimating the correlation matrix,
while reducing the confidence on the target labels with label smoothing.

##  3 PROPOSED METHOD

We construct our SER framework based on the pipeline and principles of SUPERB
[9]. As finetuning pretrained models has a high resource demand in terms of
the computational power needed, we use a simple framework with a frozen
pretrained model and lightweight classification heads. An overview of our
framework is presented in Fig. 1. In this section we describe details of the
proposed approach.

Fig. 1: Overview of the proposed architecture.

###  3.1 Layer-wise pooling

We extract information relevant for our downstream task as in SUPERB ([9]), by
taking representations from all transformer layers in the model and collapsing
them to one via a weighted average (see Fig. 1). There is one weight for each
layer (a total of $L+1$) and all weights are trained jointly with the
classification network. The weighted average representation is expressed as
follows

${\bf h}_{t}=\sum_{l=0}^{L}\gamma_{l}{\bf h}_{t,l},$ (1)

where the weights $\sum_{l=0}^{L}\gamma_{l}=1,\,\gamma_{l}\geq 0$ are
implemented with a learnable vector of size $L+1$, followed by the Softmax
function, and ${\bf h}_{t,l}$ is the representation of the $l$th layer at time
$t$ (${\bf h}_{t,0}$ is the output of the ConvNet).

###  3.2 Frame-wise pooling

Tasks requiring sentence-level classification typically employ a pooling
method, such as mean, max or attentive pooling. Mean pooling, which is
employed in SUPERB is defined as

${\bf r}=\bar{\bf h}=\frac{1}{T}\sum_{t=1}^{T}{\bf h}_{t},$ (2)

where $T$ is the number of acoustic features of an utterance extracted by the
ConvNet, ${\bf r}$ is the resulting pooled representation, while ${\bf h}_{t}$
are the representations at time $t$ after layer-wise pooling. Concatenating
the pooled representations with std features is in general helpful in speaker
recognition [23], and is implemented as

${\bf r}=\left[\bar{\bf h};\left(\frac{1}{T}\sum_{t=1}^{T}({\bf
h}_{t}-\bar{\bf h})^{2}\right)^{1/2}\right],$ (3)

where $[\cdot;\cdot]$ denotes vector concatenation and the exponents should be
considered as element-wise operators.

###  3.3 Correlation pooling

Correlation pooling (introduced in [21]) is an alternative pooling method
which has shown improvements in speaker recognition. The embedding dimension
of SSL models (typically $768$ or $1024$) is too high to estimate
correlations. We therefore project ${\bf h}$ onto a lower $d_{v}$-dimensional
space ${\bf v}$ via a linear layer ($d_{v}=256$). We then calculate the mean
vector ${\boldsymbol{\mu}}$ and the covariance matrix ${\bf\Sigma}$ of ${\bf
v}$ as follows

${\boldsymbol{\mu}}=\frac{1}{T}\sum_{t=1}^{T}{\bf
v}_{t},\,\,{\bf\Sigma}=\frac{1}{T}\sum_{t=1}^{T}({\bf
v}_{t}-{\boldsymbol{\mu}})({\bf v}_{t}-{\boldsymbol{\mu}})^{\prime},$ (4)

where ${\bf x}^{\prime}$ denotes the transpose of ${\bf x}$. Finally, the
correlation matrix is derived by normalizing with respect to the variances as
follows

${\bf C}={\bf\Sigma}\oslash{\bf S},$ (5)

where ${\bf S}={\bf s}{\bf s}^{\prime}+\epsilon{\bf 1}$, ${\bf
s}=\mbox{diag}({\bf\Sigma})^{1/2}$ (i.e. the vector of std), $\oslash$ denotes
element-wise division, while $\epsilon=10^{-8}$. Since ${\bf C}$ is a
symmetric matrix and its diagonal elements are equal to 1, we vectorize the
elements above the diagonal, yielding a $(d_{v}\times(d_{v}-1)/2)$-sized
vector, which we project onto a linear layer followed by the Softmax over the
emotion classes. For regularization, dropout is applied to ${\bf v}$, where
whole channels are dropped with probability $p_{d}=0.25$.

###  3.4 Attentive correlation pooling

We introduce here the attentive correlation pooling, by inserting weights in
the estimates of the statistics, i.e.

${\boldsymbol{\mu}}=\sum_{t=1}^{T}w_{t}{\bf
v}_{t},\,\,{\bf\Sigma}=\sum_{t=1}^{T}w_{t}({\bf
v}_{t}-{\boldsymbol{\mu}})({\bf v}_{t}-{\boldsymbol{\mu}})^{\prime}$ (6)

and we calculate ${\bf C}$ as in eq. (5). The weights $\\{w_{t}\\}_{t=1}^{T}$
(where $\sum_{t}w_{t}=1$ and $w_{t}\geq 0$) are estimated using a new flavor
of attention. Similarly to the single-head attention a single set of weights
is estimated. Similarly to the multi-head attention, multiple heads are
employed, however their similarities with ${\bf v}_{t}$ are aggregated prior
to the Softmax function via the log-sum-exp function, as follows

$\\{w_{t}\\}_{t=1}^{T}=\mbox{Softmax}\left(\\{a_{t}\\}_{t=1}^{T}\right),$ (7)

where

$a_{t}=\log\sum_{h}\exp\left({\bf q}_{h}^{\prime}{\bf o}_{t}+b_{h}\right),$
(8)

the heads $\\{{\bf q}_{h},b_{h}\\}_{h=1}^{H}$ are trainable
$d_{v}$-dimensional vectors and biases, ${\bf o}_{t}=\mbox{ReLU}({\bf
W}_{att}{\bf v}_{t})$ and ${\bf W}_{att}$ is a square matrix ($d_{v}\times
d_{v}$). Note that an equivalent implementation is to use as input to the
Softmax the $(H\times T)$-sized _vector_ of dot-products ${\bf
q}_{h}^{\prime}{\bf o}_{t}+b_{h}$ and sum the outputs over heads to obtain
$\\{w_{t}\\}_{t=1}^{T}$.

As we observe, the proposed attention resembles a mixture model with heads
parametrizing the mixture components. The log-sum-exp function is a soft
version of the max operator, meaning that $a_{t}$ is high when at least one of
the $H$ head-specific dot-products $\\{{\bf q}_{h}^{\prime}{\bf
o}_{t}+b_{h}\\}_{h=1}^{H}$ is high.

The rationale for proposing this kind of attention is two-hold. We desire to
keep the multi-modality of multi-head attention since a single head is too
weak to capture the phonetic, speaker, emotion and channel variability. On the
other hand, the standard multi-head attention results in $H$ context-vectors
(in our case $H$ correlation matrices), which can be hard to estimate
robustly, especially when the utterances are short and the estimation involves
second-order statistics.

###  3.5 Label smoothing

With label smoothing we soften the hard (one-hot) targets vectors ${\bf y}$ of
the training set as follows

${\bf y}^{LS}={\bf y}(1-p_{l})+p_{l}/K,$ (9)

where $p_{l}$ is the label smoothing parameter (i.e. the probability mass
equally distributed to all classes) and ${\bf y}^{LS}$ is the smoothed target
vector [24]. Cross-entropy is still employed as loss function, but with soft
targets.

##  4 EXPERIMENTS

###  4.1 Datasets

The IEMOCAP database consists of multi-modal recordings (speech, video) by 10
actors in dyadic sessions in English ($\approx 12$ hours) [18]. The dataset is
split in 5 dialogue sessions (one female-male speaker pair per session). The
emotions conveyed are happiness, anger, excitement, sadness, surprise, fear,
frustration, and neutral state. As in other studies on IEMOCAP, we relabel
excitement as happiness and use 4 balanced emotion classes, namely: anger,
happiness, sadness, and neutral [9, 5, 25]. All other classes are discarded.

###  4.2 Experimental Setup

We use a 5-fold cross-validation setup where at each fold we leave out one
session from the dataset. Each held-out session consists of two speakers that
are not present in the train and validation sets. This approach leaves
approximately 19% of the data for testing. Mean and standard deviation across
folds is computed and presented as the aggregated result. Our SER framework is
evaluated with WavLM, Wav2vec 2.0, and HuBERT speech representations.

##  5 RESULTS AND ANALYSIS

An overview of the results for the most common pooling methods and our
proposed approach is shown in Table 1. The results are presented for the best
configuration of our framework with $p_{d}$ and $p_{l}$ both equal to 0.25 and
$H=4$. The best overall performance was achieved for our proposed approach
using WavLM (75.60%; see also Fig. 2) which is higher compared to other SSL
approaches in the literature on the same data — 67.20% [5], 67.62% [9], 73.01%
with fine-tuned model [15].

Table 1: Unweighted accuracy (% mean and std) between test sets for SER in IEMOCAP using HuBERT large, Wav2vec 2.0, and WavLM large self-supervised representations. Pooling method | HuBERT | Wav2vec 2.0 | WavLM
---|---|---|---
mean | 65.73 (2.73) | 66.86 (1.76) | 69.44 (1.53)
mean-std | 69.15 (1.61) | 69.92 (1.17) | 72.56 (1.67)
corr ($p_{d}=0$) | 69.82 (1.35) | 68.44 (1.85) | 72.34 (1.54)
corr ($p_{d}=0.25$) | 69.72 (1.19) | 67.85 (1.84) | 72.27 (1.45)
corr attentive | 73.86 (2.10) | 70.01 (2.20) | 75.60 (2.33)

Fig. 2: Confusion matrix for WavLM using attentive correlation (left) and
meanstd pooling (right). Fig. 3: Dropout-label smoothing interaction for
Session 1 - HuBERT. Lines represent different label smoothing values.

###  5.1 Attention

To investigate the performance of the proposed attentive correlation pooling
method we experimented with different numbers of attention heads ($H$).
Specifically, we tested for $H=1,4,16,32$ heads. We obtained the best result
with $H=4$. Note that we did not observe any correlation between heads and
emotion classes, meaning that there is no direct mapping between heads and
emotions. In particular, for Session 1, HuBERT, $p_{d}=0.25,p_{l}=0.25$
accuracy was 70.32% ($H=1$), 73.18% ($H=4$), 72.53% ($H=16$), and 71.34%
($H=32$).

###  5.2 Label smoothing and dropout

To better understand the performance of our system, we probed our setup by
varying the parameters for label smoothing and dropout rate. For a set number
of attention heads ($H=4$), label smoothing was varied for
$p_{l}=0,0.15,0.25,0.3$ and dropout rate for $p_{d}=0...0.5$. The effect was
evaluated on Session 1 of the setup for attentive correlation pooling using
HuBERT and the results can be seen in Fig. 3. Both dropout and label smoothing
have an impact on the performance with label smoothing having a greater
positive effect; increasing performance from $67.93\%$ ($p_{d}=0,p_{l}=0$) to
$70,41\%$ ($p_{d}=0,p_{l}=0.25$) and even further with increasing dropout to
$73.18\%$ ($p_{d}=0.25,p_{l}=0.25$).

The impact of dropout rate and label smoothing was also investigated for mean,
mean-std, and correlation pooling. The impact in all was small to negligible.
For example, for $p_{d}=0,p_{l}=0.25$ the performance remained unchanged
compared to $p_{d}=0,p_{l}=0$ for mean, mean-std, and correlation pooling
while for $p_{d}=0.25,p_{l}=0.25$ there was a small improvement for
correlation pooling (from $68.20\%$ to $70.60\%$).

##  6 CONCLUSIONS

In this work we presented an SER framework that uses self-supervised
representations and is based on label smoothing and a novel approach to
attention, attentive correlation pooling. Notably, our method does not require
fine-tuning of the pre-trained SSL models but rather uses a light-weight
classification head that attempts to capture all relevant emotion information
from the pre-trained representations. We run several experiments using a
5-fold cross-validation setup and we have clearly demonstrated that our method
reaches high performance in all pre-trained models tested surpassing that of
the literature in similar tasks. In future work, we will extend the evaluation
setup and validate the performance of our method on more datasets.

##  7 ACKNOWLEDGEMENTS

This work was supported by the Academy of Finland project no. 340125
“Computational Modeling of Prosody in Speech”, Czech National Science
Foundation (GACR) project NEUREM3 No. 19-26934X, Czech Ministry of Interior
project No. VJ01010108 ”ROZKAZ” and Horizon 2020 Marie Sklodowska-Curie grant
ESPERANTO, No. 101007666. The authors wish to acknowledge CSC – IT Center for
Science, Finland, for providing the computational resources.

## References

  * [1] Scott Brave and Cliff Nass,  “Emotion in human-computer interaction,”  in The human-computer interaction handbook, pp. 103–118. CRC Press, 2007.
  * [2] Chul Min Lee and Shrikanth S Narayanan,  “Toward detecting emotions in spoken dialogs,”  IEEE transactions on speech and audio processing, vol. 13, no. 2, pp. 293–303, 2005.
  * [3] Rosalind W Picard,  Affective computing,  MIT press, 2000.
  * [4] Dimitrios Ververidis and Constantine Kotropoulos,  “Emotional speech recognition: Resources, features, and methods,”  Speech communication, vol. 48, no. 9, pp. 1162–1181, 2006.
  * [5] Leonardo Pepino, Pablo Riera, and Luciana Ferrer,  “Emotion recognition from speech using wav2vec 2.0 embeddings,”  arXiv preprint arXiv:2104.03502, 2021.
  * [6] Wei-Ning Hsu, Benjamin Bolte, et al.,  “HuBERT: Self-supervised speech representation learning by masked prediction of hidden units,”  IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 29, pp. 3451–3460, 2021.
  * [7] Alexei Baevski, Yuhao Zhou, et al.,  “wav2vec 2.0: A framework for self-supervised learning of speech representations,”  Advances in Neural Information Processing Systems, vol. 33, pp. 12449–12460, 2020.
  * [8] Sanyuan Chen, Chengyi Wang, et al.,  “WavLM: Large-scale self-supervised pre-training for full stack speech processing,”  arXiv preprint arXiv:2110.13900, 2021.
  * [9] Shu-wen Yang, Po-Han Chi, et al.,  “SUPERB: Speech processing universal performance benchmark,”  in Proceedings of Interspeech, 2021.
  * [10] Mousmita Sarma, Pegah Ghahremani, et al.,  “Emotion identification from raw speech signals using dnns.,”  in Interspeech, 2018, pp. 3097–3101.
  * [11] Yang Liu, Haoqin Sun, et al.,  “Discriminative feature representation based on cascaded attention network with adversarial joint loss for speech emotion recognition,”  Proc. Interspeech 2022, pp. 4750–4754, 2022.
  * [12] Jinkyu Lee and Ivan Tashev,  “High-level feature representation using recurrent neural network for speech emotion recognition,”  in Interspeech 2015, 2015.
  * [13] Yelin Kim and Emily Mower Provost,  “Leveraging inter-rater agreement for audio-visual emotion recognition,”  in 2015 International Conference on Affective Computing and Intelligent Interaction (ACII). IEEE, 2015, pp. 553–559.
  * [14] Emily Mower, Angeliki Metallinou, et al.,  “Interpreting ambiguous emotional expressions,”  in 2009 3rd International Conference on Affective Computing and Intelligent Interaction and Workshops. IEEE, 2009, pp. 1–8.
  * [15] Jiawang Liu and Haoxiang Wang,  “A speech emotion recognition framework for better discrimination of confusions.,”  in Interspeech, 2021, pp. 4483–4487.
  * [16] Haytham M Fayek, Margaret Lech, and Lawrence Cavedon,  “Modeling subjectiveness in emotion recognition with deep neural networks: Ensembles vs soft labels,”  in 2016 international joint conference on neural networks (IJCNN). IEEE, 2016, pp. 566–570.
  * [17] Lorenzo Tarantino, Philip N Garner, et al.,  “Self-attention for speech emotion recognition.,”  in Interspeech, 2019, pp. 2578–2582.
  * [18] Carlos Busso, Murtaza Bulut, et al.,  “IEMOCAP: Interactive emotional dyadic motion capture database,”  Language resources and evaluation, vol. 42, no. 4, pp. 335–359, 2008.
  * [19] Leon A Gatys, Alexander S Ecker, and Matthias Bethge,  “A neural algorithm of artistic style,”  arXiv preprint arXiv:1508.06576, 2015.
  * [20] Jan Chorowski, Ron J Weiss, et al.,  “On using backpropagation for speech texture generation and voice conversion,”  in 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2018, pp. 2256–2260.
  * [21] Themos Stafylakis, Johan Rohdin, and Lukas Burget,  “Speaker embeddings by modeling channel-wise correlations,”  in Interspeech, 2021.
  * [22] Themos Stafylakis, Ladislav Mošner, et al.,  “Extracting speaker and emotion information from self-supervised speech models via channel-wise correlations,”  in 2023 IEEE Workshop on Spoken Language Technology (SLT), to appear, pp. 1–8.
  * [23] Shuai Wang, Yexin Yang, et al.,  “Revisiting the statistics pooling layer in deep speaker embedding learning,”  in 2021 12th International Symposium on Chinese Spoken Language Processing (ISCSLP), 2021, pp. 1–5.
  * [24] Christian Szegedy, Vincent Vanhoucke, et al.,  “Rethinking the inception architecture for computer vision,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 2818–2826.
  * [25] Haytham M Fayek, Margaret Lech, and Lawrence Cavedon,  “Evaluating deep learning architectures for speech emotion recognition,”  Neural Networks, vol. 92, pp. 60–68, 2017.


11institutetext:  School of Basic Sciences (Mathematics), Indian Institute of
Technology Bhubaneswar, India. 11email<EMAIL_ADDRESS>
# Convergence of the Neumann-Neumann Method for the Cahn-Hilliard Equation

Gobinda Garai 11

###### Abstract

In this paper, we analyze a substructuring type algorithm for the Cahn-
Hilliard (CH) equation. Being a nonlinear equation, it is of great importance
to develop efficient numerical schemes for investigating the solution
behaviour of the CH equation. We present the formulation of Neumann-Neumann
(NN) method applied to the CH equation and investigate the convergence
behaviour of the same in one and two spatial dimension for two subdomains. We
illustrate the theoretical results by providing numerical example.

###### keywords:

Neumann-Neumann, Domain Decomposition, Parallel computing, Cahn-Hilliard
equation

##  1 Introduction

The Cahn-Hilliard equation

$\displaystyle\frac{\partial u}{\partial t}$ $\displaystyle=\Delta v,\,\
\mbox{ for}\,\ (x,t)\in\Omega\times(0,T],$ (1) $\displaystyle v$
$\displaystyle=F^{\prime}(u)-\epsilon^{2}\Delta u,\,\ \mbox{ for}\,\
(x,t)\in\Omega\times(0,T],$ $\displaystyle u(x,0)$ $\displaystyle=u_{0}(x),\,\
x\in\Omega,$

where $\Omega\subset\mathbb{R}^{d}(d=1,2)$, with the Neumann boundary
condition $\frac{\partial u}{\partial n}=\frac{\partial v}{\partial n}=0\mbox{
on}\,\ \partial\Omega,$ for all $t\in(0,T]$, has been suggested to represent
the evolution of a binary melted alloy below the critical temperature in [4,
5]. The solution of (1) involves two different dynamics, one is phase
separation which is quick in time, and another is phase coarsening which is
relatively slower in time. The fine-scale phase regions are emerged in the
early period and separated by the interface which is of width $\epsilon$.
Whereas during phase coarsening, the solution lean toward an equilibrium state
which reduces the internal energy. This means, the equation (1) describes the
total mass conservation and energy minimization while the system evolves. The
energy minimization is a key feature of (1) and is expected to be preserved by
numerical method and to deal with that, Eyre [6, 7] proposed an
unconditionally gradient stable scheme. The idea is to split the homogeneous
free energy $F(u)=0.25(u^{2}-1)^{2}$ into the sum of a convex and a concave
term, and then treating the convex term implicitly and the concave term
explicitly to obtain, for example, a first order in time and 2nd order in
space approximation for (1) in 1D:

$\displaystyle u_{j}^{n+1}$
$\displaystyle=\delta\Delta_{h}v_{j}^{n+1}+u_{j}^{n},$ (2) $\displaystyle
v_{j}^{n+1}$
$\displaystyle=(u_{j}^{n+1})^{3}-u_{j}^{n}-\epsilon^{2}\Delta_{h}u_{j}^{n+1},$

where $\delta$ is the time step and $\Delta_{h}$ is the discrete Laplacian.
The equation (2) represents a system of nonlinear coupled equation due to the
cubic term. To linearise the problem, the term $(u_{j}^{n+1})^{3}$ is
rewritten as $(u_{j}^{n})^{2}u_{j}^{n+1}$, which results the following
modified system

$\displaystyle u_{j}^{n+1}-\delta\Delta_{h}v_{j}^{n+1}$
$\displaystyle=u_{j}^{n},$ $\displaystyle
v_{j}^{n+1}+\epsilon^{2}\Delta_{h}u_{j}^{n+1}-(u_{j}^{n})^{2}u_{j}^{n+1}$
$\displaystyle=-u_{j}^{n},$

which is also an unconditionally gradient stable scheme and has the same
accuracy as the nonlinear scheme (2) [6]. Therefore one has to solve the
following system of elliptic equations at each time level

$\displaystyle\begin{bmatrix}I&-\delta\Delta\\\
\epsilon^{2}\Delta-c^{2}&I\end{bmatrix}\quad\begin{bmatrix}{\bar{u}}\\\
{\bar{v}}\end{bmatrix}\quad=\quad\begin{bmatrix}f_{\bar{u}}\\\
f_{\bar{v}}\end{bmatrix},\quad\mbox{in}\,\ \Omega,$ (3)

where $c=u_{j}^{n},\bar{u}=u_{j}^{n+1},\bar{v}=v_{j}^{n+1},f_{\bar{u}}=c\,\
\text{and}\,\ f_{\bar{v}}=-c.$ In 2D, one also gets the above system (3) at
each time level for suitably chosen $f_{\bar{u}},f_{\bar{v}},c$. It is worth
mentioning that many other basic algorithms approximating the solution of the
CH equation can be reformulated as (3), for example the semi-implicit Euler’s
scheme [9], the Linearly stabilized splitting scheme [1].

Since the spatial mesh size $h$ is $O(\epsilon)$ or even finer, the linear
system (3) will result in a massive algebraic structure and which needs to be
solved serially for getting the long term solution of (1). Thus, it is of
great importance to fast-track the simulation using parallel computation,
which can be fulfilled by domain decomposition techniques [10, 11]. In this
work, we lay our efforts on the Neumann-Neumann method [3, 8, 2]. The main
objective of our work is to solve the problem (3) with the imposed
transmission condition and analyze the convergence behaviour for two subdomain
setting in 1D and 2D.

We introduce the NN algorithm in one and two spatial dimension for two
subdomains, and state our convergence result in Section 2. To illustrate our
analysis, the accuracy and robustness of the proposed formulation, we present
the numerical results in Section 3.

##  2 The Neumann Neumann Method

In this section, we introduce the Neumann-Neumann method for the second order
elliptic system (3). For convenience we use the notation $u,v$ instead of
$\bar{u},\bar{v}$ and rewrite the system (3) as

$\displaystyle\begin{bmatrix}I&-\delta\Delta\\\
\epsilon^{2}\Delta-c^{2}&I\end{bmatrix}\begin{bmatrix}{u}\\\
{v}\end{bmatrix}=\begin{bmatrix}f_{{u}}\\\
f_{{v}}\end{bmatrix},\quad\mbox{in}\,\ \Omega,$ (4)

together with the physical boundary condition $\mathcal{B}\begin{bmatrix}u\\\
v\end{bmatrix}=0.$

We present the NN method for the problem (4) on the spatial domain $\Omega$
with Neumann boundary conditions on $\partial\Omega$. Suppose the spatial
domain $\Omega$ is partitioned into two non-overlapping subdomains
$\Omega_{1},\Omega_{2}$. We denote $u_{j}$ the restriction of the solution $u$
of (4) to $\Omega_{j}$ for $j=1,2$ and set
$\Gamma:=\partial\Omega_{1}\cap\partial\Omega_{2}$.

The NN algorithm for the CH system (4) starts with initial guesses
$g^{[0]},h^{[0]}$ along the interface $\Gamma$ and solve for $k=1,2,...$

$\displaystyle\left\\{\begin{aligned} \begin{bmatrix}I&-\delta\Delta\\\
\epsilon^{2}\Delta-c^{2}&I\end{bmatrix}\begin{bmatrix}u_{j}^{[k]}\\\
v_{j}^{[k]}\end{bmatrix}=\begin{bmatrix}f_{u}\\\ f_{v}\end{bmatrix},\,\
\text{in}\,\ \Omega_{j},\\\ \mathcal{B}\begin{bmatrix}u_{j}^{[k]}\\\
v_{j}^{[k]}\end{bmatrix}=0,\,\ \mbox{on}\,\
\partial\Omega_{j}\cap\partial\Omega,\\\ \begin{bmatrix}u_{j}^{[k]}\\\
v_{j}^{[k]}\end{bmatrix}=\begin{bmatrix}g^{[k-1]}\\\
h^{[k-1]}\end{bmatrix},\,\ \mbox{on}\,\ \Gamma,\\\ \end{aligned}\right.$
$\displaystyle\left\\{\begin{aligned} \begin{bmatrix}I&-\delta\Delta\\\
\epsilon^{2}\Delta-c^{2}&I\end{bmatrix}\begin{bmatrix}\phi_{j}^{[k]}\\\
\psi_{j}^{[k]}\end{bmatrix}=0,\,\ \text{in}\,\ \Omega_{j},\\\
\mathcal{B}\begin{bmatrix}\phi_{j}^{[k]}\\\ \psi_{j}^{[k]}\end{bmatrix}=0,\,\
\mbox{on}\,\ \partial\Omega_{j}\cap\partial\Omega,\\\ \frac{\partial}{\partial
n}\begin{bmatrix}\phi_{j}^{[k]}\\\
\psi_{j}^{[k]}\end{bmatrix}=\frac{\partial}{\partial
n}\begin{bmatrix}u_{1}^{[k]}-u_{2}^{[k]}\\\
v_{1}^{[k]}-v_{2}^{[k]}\end{bmatrix},\,\ \mbox{on}\,\ \Gamma.\\\
\end{aligned}\right.$ (5)

Then we update the interface trace by

$\begin{bmatrix}g^{[k]}\\\ h^{[k]}\end{bmatrix}=\begin{bmatrix}g^{[k-1]}\\\
h^{[k-1]}\end{bmatrix}-\theta\begin{bmatrix}\phi_{1}^{[k]}-\phi_{2}^{[k]}\\\
\psi_{1}^{[k]}-\psi_{2}^{[k]}\end{bmatrix}_{\big{|}_{\Gamma}},$

where $\theta\in(0,1)$ is the relaxation parameter. For the convergence
analysis we consider the error equation, which is a homogeneous version of (5)
by the virtue of linearity. The goal of our study is to analyze how the error
$g^{[k]},h^{[k]}$ converges to zero as $k\rightarrow\infty$.

###  2.1 Convergence analysis in 1D

Let the spatial domain $\Omega=(-a,b)$ is decomposed into $\Omega_{1}=(-a,0)$
and $\Omega_{2}=(0,b)$ with interface $\Gamma=\\{0\\}$. We consider the error
equation described in (5) by taking $f_{u}=0,f_{v}=0$. The system (5) can be
rewritten in 1D as:

$\displaystyle\left\\{\begin{aligned} AE_{j}^{[k]}&=0,\quad\quad\text{in}\,\
\Omega_{j}\\\ E_{j}^{[k]}&=\begin{bmatrix}g^{[k-1]}\\\
h^{[k-1]}\end{bmatrix},\quad\text{on}\,\ \Gamma\\\ \frac{\partial}{\partial
x}\ E_{j}^{[k]}&=0,\quad\quad\text{on}\,\
\partial\Omega_{j}\cap\partial\Omega\\\ \end{aligned}\right.$
$\displaystyle\left\\{\begin{aligned} AF_{j}^{[k]}&=0,\quad\quad\text{in}\,\
\Omega_{i}\\\ \frac{\partial}{\partial x}\
F_{j}^{[k]}&=\frac{\partial}{\partial
x}\left[E_{1}^{[k]}-E_{2}^{[k]}\right],\quad\text{on}\,\ \Gamma\\\
\frac{\partial}{\partial x}\ F_{j}^{[k]}&=0,\quad\text{on}\,\
\partial\Omega_{j}\cap\partial\Omega\\\ \end{aligned}\right.$ (6)

with the update condition

$\begin{bmatrix}g^{[k]}\\\ h^{[k]}\end{bmatrix}=\begin{bmatrix}g^{[k-1]}\\\
h^{[k-1]}\end{bmatrix}-\theta\left[F_{1}^{[k]}-F_{2}^{[k]}\right]_{\big{|}_{\Gamma}},$
(7)

where

$A=\begin{bmatrix}1&-\delta\frac{d^{2}}{dx^{2}}\\\
\epsilon^{2}\frac{d^{2}}{dx^{2}}-c^{2}&1\end{bmatrix},\quad
E_{j}^{[k]}=\begin{bmatrix}u_{j}^{[k]}(x)\\\ v_{j}^{[k]}(x)\end{bmatrix},\quad
F_{j}^{[k]}=\begin{bmatrix}\phi_{j}^{[k]}(x)\\\
\psi_{j}^{[k]}(x)\end{bmatrix},$

where $E_{j},F_{j}$ are the subdomain solutions in $\Omega_{j}$ for $j=1,2$ at
the $k-$th iteration.

To analyze the convergence of the algorithm (5), we need to determine the
subdomain solutions by solving the following algebraic equations

$AE_{j}=0,j=1,2.$ (8)

We assume the solution of the equation (8) is of the following form

$E_{j}=\Psi_{j}e^{\xi x},\text{where}\;\xi\;\text{is parameter to be
determined.}$ (9)

Inserting the above expression of $E_{j}$ into equation (8) and using the non-
vanishing property of exponential term, we get

$\begin{bmatrix}1&-\delta\xi^{2}\\\
\epsilon^{2}\xi^{2}-c^{2}&1\end{bmatrix}\Psi_{j}=0.$ (10)

Now the equation (10) has non-trivial solutions only if the coefficient matrix
is singular, i.e., the determinant of the coefficient matrix is zero,

$\det\begin{bmatrix}1&-\delta\xi^{2}\\\
\epsilon^{2}\xi^{2}-c^{2}&1\end{bmatrix}=0.$

Solving this equation yields
$\xi_{1,2}=\pm\sqrt{\lambda_{1}},\quad\xi_{3,4}=\pm\sqrt{\lambda_{2}}$, where
$\lambda_{1,2}$ are given by

$\lambda_{1,2}=\frac{c^{2}\delta\pm\sqrt{c^{4}\delta^{2}-4\epsilon^{2}\delta}}{2\epsilon^{2}\delta}$

respectively. Thus the solution of (8) has the following form

$E_{j}^{[k]}=\sum_{l=1}^{4}\zeta_{j,l}^{[k]}\mu_{l}e^{\xi_{l}x},$

where for each $l$, $\zeta_{j,l}$ are constant weights for $j=1,2$ and
$\mu_{l}$ is the eigenvector of the coefficient matrix in (10) associated to
the eigenvalue zero for $\xi=\xi_{l}$, and is explicitly given by

$\mu_{1}=\mu_{2}=\begin{bmatrix}\delta\lambda_{1}\\\
1\end{bmatrix},\quad\mu_{3}=\mu_{4}=\begin{bmatrix}\delta\lambda_{2}\\\
1\end{bmatrix}.$

Using the transmission conditions on the interface $\Gamma$ and the physical
boundary conditions on $\partial\Omega_{i}\backslash\Gamma$, we can get rid of
the constants $\zeta_{j,l}$ for each subdomain. Set
$\sigma_{i,x}:=\sinh(\xi_{i}x),\sigma_{i,a}:=\sinh(\xi_{i}a),\sigma_{i,b}:=\sinh(\xi_{i}b),\gamma_{i,x}:=\cosh(\xi_{i}x),\gamma_{i,a}:=\cosh(\xi_{i}a),\gamma_{i,b}:=\cosh(\xi_{i}b)$
for $i=1,3$ and
$\eta_{1}=\frac{g^{k-1}-\delta\lambda_{2}h^{k-1}}{\delta\lambda},\eta_{2}=\frac{g^{k-1}-\delta\lambda_{1}h^{k-1}}{-\delta\lambda}$
with $\lambda=\lambda_{1}-\lambda_{2}$, and
$\nu_{i}=\frac{\sigma_{i,a}}{\gamma_{i,a}}+\frac{\sigma_{i,b}}{\gamma_{i,b}}$
for $i=1,3$. We have the subdomain solution at $k-$th iteration for Dirichlet
step

$E_{1}^{[k]}=\begin{bmatrix}\mu_{1}&\mu_{3}\end{bmatrix}\begin{bmatrix}\gamma_{1,x}+\frac{\sigma_{1,a}\sigma_{1,x}}{\gamma_{1,a}}&0\\\
0&\gamma_{3,x}+\frac{\sigma_{3,a}\sigma_{3,x}}{\gamma_{3,a}}\end{bmatrix}\begin{bmatrix}\eta_{1}\\\
\eta_{2}\end{bmatrix},$

$E_{2}^{[k]}=\begin{bmatrix}\mu_{1}&\mu_{3}\end{bmatrix}\begin{bmatrix}\gamma_{1,x}-\frac{\sigma_{1,b}\sigma_{1,x}}{\gamma_{1,b}}&0\\\
0&\gamma_{3,x}-\frac{\sigma_{3,b}\sigma_{3,x}}{\gamma_{3,b}}\end{bmatrix}\begin{bmatrix}\eta_{1}\\\
\eta_{2}\end{bmatrix},$

and for Neumann step we solve $AF_{j}=0$ for $j=1,2$ with physical boundary
conditions on $\partial\Omega_{i}\backslash\Gamma$ and imposed transmission
condition on the interface $\Gamma$ and get the subdomain solution at $k-$th
iteration,

$F_{1}^{[k]}=\begin{bmatrix}\mu_{1}&\mu_{3}\end{bmatrix}\begin{bmatrix}\sigma_{1,x}\nu_{1}+\frac{\gamma_{1,x}\gamma_{1,a}\nu_{1}}{\sigma_{1,a}}&0\\\
0&\sigma_{3,x}\nu_{3}+\frac{\gamma_{3,x}\gamma_{3,a}\nu_{3}}{\sigma_{3,a}}\end{bmatrix}\begin{bmatrix}\eta_{1}\\\
\eta_{2}\end{bmatrix},$

$F_{2}^{[k]}=\begin{bmatrix}\mu_{1}&\mu_{3}\end{bmatrix}\begin{bmatrix}\sigma_{1,x}\nu_{1}-\frac{\gamma_{1,x}\gamma_{1,b}\nu_{1}}{\sigma_{1,b}}&0\\\
0&\sigma_{3,x}\nu_{3}-\frac{\gamma_{3,x}\gamma_{3,b}\nu_{3}}{\sigma_{3,b}}\end{bmatrix}\begin{bmatrix}\eta_{1}\\\
\eta_{2}\end{bmatrix}.$

Using the subdomain solutions in (7) we get the recurrence relation

$\begin{bmatrix}g^{[k]}\\\
h^{[k]}\end{bmatrix}=\mathbb{P}\begin{bmatrix}g^{[k-1]}\\\
h^{[k-1]}\end{bmatrix},$ (11)

where the iteration matrix $\mathbb{P}$ is given by

$\mathbb{P}=\begin{bmatrix}1-\theta\frac{\lambda_{1}\rho_{1}-\lambda_{2}\rho_{2}}{\lambda}&\theta\frac{\delta\lambda_{1}\lambda_{2}(\rho_{1}-\rho_{2})}{\lambda}\\\
\\\
\theta\frac{\rho_{2}-\rho_{1}}{\delta\lambda}&1-\theta\frac{\lambda_{1}\rho_{2}-\lambda_{2}\rho_{1}}{\lambda}\end{bmatrix},$

with
$\rho_{1}=2+\frac{\sigma_{1,a}\gamma_{1,b}}{\gamma_{1,a}\sigma_{1,b}}+\frac{\sigma_{1,b}\gamma_{1,a}}{\gamma_{1,b}\sigma_{1,a}},\;\rho_{2}=2+\frac{\sigma_{3,a}\gamma_{3,b}}{\gamma_{3,a}\sigma_{3,b}}+\frac{\sigma_{3,b}\gamma_{3,a}}{\gamma_{3,b}\sigma_{3,b}}$.

######  Theorem 2.1 (Convergence of NN for symmetric case).

When the subdomains are of same size, $a=b$, the NN algorithm for two
subdomains converges linearly for $0<\theta<1/2,\theta\neq\ 1/4$. Moreover for
$\theta=1/4$, it converges in two iterations.

######  Proof 2.2.

If $a=b$, then (11) becomes

$\begin{bmatrix}g^{[k]}\\\ h^{[k]}\end{bmatrix}=\begin{bmatrix}1-4\theta&0\\\
0&1-4\theta\end{bmatrix}\begin{bmatrix}g^{[k-1]}\\\ h^{[k-1]}\end{bmatrix},$

which means the method converges to the exact solution in two iterations for
$\theta=1/4$. Also, the algorithm diverges if $|1-4\theta|\geq 1$, i.e if
$\theta\geq 1/2$. Thus the convergence is linear for $0<\theta<1/2,\theta\neq
1/4$.

We now study the convergence results for unequal sudomains with the special
value of $\theta$ being $1/4$. Before going to the main result, we prove the
following Lemma, which is needed for future analysis.

######  Lemma 2.3.

The function $f(t)=\frac{\sinh^{2}((a-b)t)}{\sinh(2at)\sinh(2bt)}$ with
$a,b>0$ has the following properties

  1. 1.

$\forall t>0,0<f(t)<\frac{(a-b)^{2}}{4ab}$.

  2. 2.

$\forall t>0,f(t)$ is a monotonically decreasing function.

######  Proof 2.4.

Clearly $f(t)>0,\forall t>0$. We prove (2) and the upper bound estimate in (1)
will follow from that. Taking the derivative of logarithm of $f(t)$ with
respect to $t$ we get

$\frac{d}{dt}(\ln(f(t)))=2(a-b)\coth((a-b)t)-2a\coth(2at)-2b\coth(2bt).$

Now $f^{\prime}(t)<0$ iff $(a-b)\coth((a-b)t)<a\coth(2at)+b\coth(2bt)$, which
is true, since the function $t\coth(t)$ is strictly convex and increasing for
$t>0$. This completes the result.

######  Theorem 2.5 (Convergence of NN for non-symmetric case).

For $\theta=1/4$ the error of the NN method for two subdomains satisfies,

$\parallel
g^{[k]}\parallel_{L^{\infty}(\Gamma)}\leq\begin{cases}\left(\frac{(a-b)^{2}}{4ab}\right)^{k}\max\left\\{\parallel
g^{[0]}\parallel_{L^{\infty}(\Gamma)},\parallel
h^{[0]}\parallel_{L^{\infty}(\Gamma)}\right\\},&\text{if}\;\delta>\frac{4\epsilon^{2}}{c^{4}},\\\
\left(\frac{(a-b)^{2}}{2ab}\right)^{k}\max\left\\{\parallel
g^{[0]}\parallel_{L^{\infty}(\Gamma)},\parallel
h^{[0]}\parallel_{L^{\infty}(\Gamma)}\right\\},&\text{if}\;\delta<\frac{4\epsilon^{2}}{c^{4}}.\end{cases}$

######  Proof 2.6.

When $\delta>\frac{4\epsilon^{2}}{c^{4}}$, it is clear that $\lambda_{1,2}$
are real and positive, so are $\xi_{1,3}$. The iteration matrix $\mathbb{P}$
in (11) has eigenvalues $(1-\frac{\rho_{1}}{4}),(1-\frac{\rho_{2}}{4})$.
Hence, the convergence is reached iff the spectral radius
$\rho(\mathbb{P})=\max\\{|1-\frac{\rho_{1}}{4}|,|1-\frac{\rho_{2}}{4}|\\}$ is
less than one. Upon simplification, we get the expression of the eigenvalues
as: $-\frac{\sinh^{2}((a-b)\xi_{i})}{\sinh(2a\xi_{i})\sinh(2b\xi_{i})}$, for
$i=1,3$. Using the Lemma (2.3) we have the required estimate.
For $\delta<\frac{4\epsilon^{2}}{c^{4}}$, $\lambda_{1}$ and $\lambda_{2}$ are
complex conjugate numbers. We denote $\lambda_{1,2}$ as
$\lambda_{1,2}=\lambda_{\Re}\pm i\lambda_{\Im},i=\sqrt{-1}$, where

$\lambda_{\Re}=\frac{c^{2}}{2\epsilon^{2}},\;\lambda_{\Im}=\frac{\sqrt{4\delta\epsilon^{2}-\delta^{2}c^{4}}}{2\delta\epsilon^{2}}.$

Then $\xi_{1,3}$ takes the form
$\xi_{1}=\xi_{\Re}+i\xi_{\Im},\xi_{3}=\xi_{\Re}-i\xi_{\Im}$, where

$\xi_{\Re}=\sqrt{\frac{\lambda_{\Re}+\sqrt{\lambda_{\Re}^{2}+\lambda_{\Im}^{2}}}{2}},\;\xi_{\Im}=\sqrt{\frac{-\lambda_{\Re}+\sqrt{\lambda_{\Re}^{2}+\lambda_{\Im}^{2}}}{2}}.$

Clearly, $\xi_{\Re},\xi_{\Im}$ are positive number. Now if we take the modulus
of the complex eigenvalues of the iteration matrix $\mathbb{P}$, we have

$\displaystyle\left|-\frac{\sinh^{2}((a-b)\xi_{i})}{\sinh(2a\xi_{i})\sinh(2b\xi_{i})}\right|$
$\displaystyle\leq$
$\displaystyle\frac{\sinh^{2}((a-b)\xi_{\Re})+\sin^{2}((a-b)\xi_{\Im})}{\sinh(2a\xi_{\Re})\sinh(2b\xi_{\Re})}$
$\displaystyle\leq$
$\displaystyle\frac{2\sinh^{2}((a-b)\xi_{\Re})}{\sinh(2a\xi_{\Re})\sinh(2b\xi_{\Re})}\leq\left(\frac{(a-b)^{2}}{2ab}\right),$

the second inequality follows from
$\sin^{2}((a-b)\xi_{\Im})<\sinh^{2}((a-b)\xi_{\Im})<\sinh^{2}((a-b)\xi_{\Re})$,
as $\xi_{\Im}<\xi_{\Re}$, and the last inequality follows from the Lemma
(2.3). Hence the estimate.

######  Remark 2.7.

The above linear estimate will imply convergence if $\frac{(a-b)^{2}}{4ab}<1$,
i.e., $q^{2}-6q+1<0$, where $q=a/b$. Hence $q$ must belongs to
$\left(\frac{6-\sqrt{32}}{2},\frac{6+\sqrt{32}}{2}\right)$ when
$\delta>\frac{4\epsilon^{2}}{c^{4}}$ and for
$\delta<\frac{4\epsilon^{2}}{c^{4}}$, $q$ is in $(2-\sqrt{3},2+\sqrt{3})$.

###  2.2 Convergence analysis in 2D

We now analyze the convergence of the NN method for a decomposition with two
subdomains in 2D. Let the domain $\Omega=(-a,b)\times(0,L)$ be partitioned
into two rectangles, given by $\Omega_{1}=(-a,0)\times(0,L)$ and
$\Omega_{2}=(0,b)\times(0,L)$. We analyze this 2D case by transforming it into
a collection of 1D problems using the Fourier sine transform along
$y$-direction. Expanding the solution
$u_{i}^{[k]},v_{i}^{[k]},\phi_{i}^{[k]},\psi_{i}^{[k]}$ for $i=1,2$ in a
Fourier sine series along $y$-direction yields

$\begin{array}[]{cc}u_{i}^{[k]}(x,y)=\sum_{m\geq
1}\hat{u}_{i}^{[k]}(x,m)\sin(\frac{m\pi y}{L}),&v_{i}^{[k]}(x,y)=\sum_{m\geq
1}\hat{v}_{i}^{[k]}(x,m)\sin(\frac{m\pi y}{L}),\\\
\phi_{i}^{[k]}(x,y)=\sum_{m\geq 1}\hat{\phi}_{i}^{[k]}(x,m)\sin(\frac{m\pi
y}{L}),&\psi_{i}^{[k]}(x,y)=\sum_{m\geq
1}\hat{\psi}_{i}^{[k]}(x,m)\sin(\frac{m\pi y}{L}).\end{array}$

After a Fourier sine transform, the NN algorithm (5) for the error equations
for 2D CH equation becomes

$\displaystyle\left\\{\begin{aligned}
\hat{A}\hat{E}_{1}^{[k]}&=0,\quad\quad\text{in}\,\ (-a,0),\\\
\hat{E}_{1}^{[k]}&=\begin{bmatrix}\hat{g}^{[k-1]}\\\
\hat{h}^{[k-1]}\end{bmatrix},\quad\text{at}\,\ x=0,\\\
\frac{\partial}{\partial x}\ \hat{E}_{1}^{[k]}&=0,\quad\quad\text{at}\,\
x=-a,\\\ \end{aligned}\right.$ $\displaystyle\left\\{\begin{aligned}
\hat{A}E_{2}^{[k]}&=0,\quad\quad\text{in}\,\ (0,b),\\\
\hat{E}_{2}^{[k]}&=\begin{bmatrix}\hat{g}^{[k-1]}\\\
\hat{h}^{[k-1]}\end{bmatrix},\quad\text{at}\,\ x=0,\\\
\frac{\partial}{\partial x}\ \hat{E}_{2}^{[k]}&=0,\quad\quad\text{at}\,\
x=b,\\\ \end{aligned}\right.$ (12)

and

$\displaystyle\left\\{\begin{aligned}
\hat{A}\hat{F}_{1}^{[k]}&=0,\quad\quad\text{in}\,\ (-a,0),\\\
\frac{\partial}{\partial x}\ \hat{F}_{1}^{[k]}&=\frac{\partial}{\partial
x}\left[\hat{E}_{1}^{[k]}-\hat{E}_{2}^{[k]}\right],\,\text{at}\,\ x=0,\\\
\frac{\partial}{\partial x}\ \hat{F}_{1}^{[k]}&=0,\quad\quad\text{at}\,\
x=-a,\\\ \end{aligned}\right.$ $\displaystyle\left\\{\begin{aligned}
\hat{A}\hat{F}_{2}^{[k]}&=0,\quad\quad\text{in}\,\ (0,b),\\\
\frac{\partial}{\partial x}\ \hat{F}_{2}^{[k]}&=\frac{\partial}{\partial
x}\left[\hat{E}_{1}^{[k]}-\hat{E}_{2}^{[k]}\right],\,\text{at}\,\ x=0,\\\
\frac{\partial}{\partial x}\ \hat{F}_{2}^{[k]}&=0,\quad\quad\text{at}\,\
x=b,\\\ \end{aligned}\right.$ (13)

and the update condition becomes

$\begin{bmatrix}\hat{g}^{[k]}\\\
\hat{h}^{[k]}\end{bmatrix}=\begin{bmatrix}\hat{g}^{[k-1]}\\\
\hat{h}^{[k-1]}\end{bmatrix}-\theta\left[\hat{F}_{1}^{[k]}-\hat{F}_{2}^{[k]}\right]_{\big{|}\\{x=0\\}},$

where

$\hat{A}=\begin{bmatrix}1&-\delta(\frac{d^{2}}{dx^{2}}-p_{m}^{2})\\\
\epsilon^{2}(\frac{d^{2}}{dx^{2}}-p_{m}^{2})-c^{2}&1\end{bmatrix},\;\hat{E}_{j}^{[k]}=\begin{bmatrix}\hat{u}_{j}^{[k]}(x,m)\\\
\hat{v}_{j}^{[k]}(x,m)\end{bmatrix},\;\hat{F}_{j}^{[k]}=\begin{bmatrix}\hat{\phi}_{j}^{[k]}(x,m)\\\
\hat{\psi}_{j}^{[k]}(x,m)\end{bmatrix},$

where $\hat{E}_{j}^{[k]}(x,m),\hat{F}_{j}^{[k]}(x,m)$ denote the Fourier sine
coefficients of $E_{j}^{[k]}(x,y)$ and $F_{j}^{[k]}(x,y)$, and
$p_{m}^{2}=\frac{\pi^{2}m^{2}}{L^{2}}$. We perform similar calculation as in
the case of one dimensional analysis and get the recurrence relation as:

$\begin{bmatrix}\hat{g}^{[k]}\\\
\hat{h}^{[k]}\end{bmatrix}=\hat{\mathbb{P}}\begin{bmatrix}\hat{g}^{[k-1]}\\\
\hat{h}^{[k-1]}\end{bmatrix},$ (14)

where $\hat{\mathbb{P}}$ is the iteration matrix given by

$\hat{\mathbb{P}}=\begin{bmatrix}1-\theta\frac{\lambda_{1}\rho_{1}-\lambda_{2}\rho_{2}}{\lambda}&\theta\frac{\delta\lambda_{1}\lambda_{2}(\rho_{1}-\rho_{2})}{\lambda}\\\
\\\
\theta\frac{\rho_{2}-\rho_{1}}{\delta\lambda}&1-\theta\frac{\lambda_{1}\rho_{2}-\lambda_{2}\rho_{1}}{\lambda}\end{bmatrix}$

with the expressions $\rho_{1},\rho_{2}$ exactly as defined earlier, except
that $\lambda_{i}$ is replaced by $\lambda_{i}+p_{m}^{2}$ in the expression of
$\xi_{1,2}$ and $\xi_{3,4}$.

######  Theorem 2.8 (Convergence of NN in 2D).

When the subdomains are of same size, $a=b$, the NN algorithm converges
linearly for $0<\theta<1/2,\theta\neq\ 1/4$, and for $\theta=1/4$, it
converges in two iterations. When the subdomains are unequal, for $\theta=1/4$
the error of the NN method for two subdomains satisfies,

$\parallel
g^{[k]}\parallel_{L^{2}(\Gamma)}\leq\begin{cases}\left(\frac{(a-b)^{2}}{4ab}\right)^{k}\max\left\\{\parallel
g^{[0]}\parallel_{L^{2}(\Gamma)},\parallel
h^{[0]}\parallel_{L^{2}(\Gamma)}\right\\},&\text{if}\;\delta>\frac{4\epsilon^{2}}{c^{4}},\\\
\left(\frac{(a-b)^{2}}{2ab}\right)^{k}\max\left\\{\parallel
g^{[0]}\parallel_{L^{2}(\Gamma)},\parallel
h^{[0]}\parallel_{L^{2}(\Gamma)}\right\\},&\text{if}\;\delta<\frac{4\epsilon^{2}}{c^{4}}.\end{cases}$

######  Proof 2.9.

If $a=b$, then (14) yields

$\begin{bmatrix}\hat{g}^{[k]}\\\
\hat{h}^{[k]}\end{bmatrix}=\begin{bmatrix}1-4\theta&0\\\
0&1-4\theta\end{bmatrix}\begin{bmatrix}\hat{g}^{[k-1]}\\\
\hat{h}^{[k-1]}\end{bmatrix},$

which upon back Fourier transform shows that the method converges to the exact
solution in two iterations for $\theta=1/4$. Also, the algorithm diverges if
$|1-4\theta|\geq 1$, i.e if $\theta\geq 1/2$. Thus the convergence is linear
for $0<\theta<1/2,\theta\neq 1/4$.

Now in the case of unequal subdomain with
$\delta>\frac{4\epsilon^{2}}{c^{4}}$, it is clear that $\lambda_{1,2}$ are
real and positive, so are $\xi_{1,3}$. The iteration matrix $\hat{\mathbb{P}}$
in (14) has spectral radius
$\rho(\hat{\mathbb{P}})=\max\left\\{|1-\frac{\rho_{1}}{4}|,|1-\frac{\rho_{2}}{4}|\right\\}$.
By the recurrence relation (14) and using Parseval-Plancherel identity we get

$\begin{array}[]{cc}\parallel
g^{[k]}\parallel_{L^{2}(\Gamma)}\leq\rho(\hat{\mathbb{P}})\max\left\\{\parallel
g^{[k-1]}\parallel_{L^{2}(\Gamma)},\parallel
h^{[k-1]}\parallel_{L^{2}(\Gamma)}\right\\},\\\ \parallel
h^{[k]}\parallel_{L^{2}(\Gamma)}\leq\rho(\hat{\mathbb{P}})\max\left\\{\parallel
g^{[k-1]}\parallel_{L^{2}(\Gamma)},\parallel
h^{[k-1]}\parallel_{L^{2}(\Gamma)}\right\\}.\end{array}$

We get the required result by estimating $\rho(\hat{\mathbb{P}})$ similar to
Theorem (2.5).
The case $\delta<\frac{4\epsilon^{2}}{c^{4}}$ can be treated similarly as in
Theorem (2.5).

##  3 Numerical Illustration

First, we present numerical results for the NN algorithm with small time
window. The parameter $\epsilon$ is taken as $0.01$, except otherwise stated.
The iterations start from a random initial guess and stop as the error reaches
a tolerance of $10^{-6}$ in $L^{\infty}\;\text{or}\;L^{2}$ norm. The phase
separation process is rapid in time, and consequently small time steps should
be taken. We choose $\delta=1\times 10^{-6}$, which falls in the bracket
$\delta<\frac{4\epsilon^{2}}{c^{4}}$. The iteration count for the NN method
are shown in Table 1 for 1D. In Fig 1, on the left we plot the error curves
for NN method for different values of $\theta$ for CH equation in 2D, and on
the right we show the numerical error and theoretical estimates for
$\theta=1/4$.

$h$ $\theta$ |  0.l |  0.2 |  0.25 |  0.3 |  0.4 |  0.5 |  0.6 |  0.7 |  0.8
---|---|---|---|---|---|---|---|---|---
1/64 |  35 |  12 |  2 |  12 |  37 |  - |  - |  - |  -
1/128 |  35 |  12 |  2 |  12 |  38 |  - |  - |  - |  -
1/256 |  36 |  12 |  2 |  12 |  39 |  - |  - |  - |  -
1/512 |  37 |  12 |  2 |  12 |  39 |  - |  - |  - |  -
Table 1: Number of iteration compared for unequal subdomain in 1D.

(a)

(b)

Figure 1: Iteration compared in 2D (left), and comparison of theoretical and
numerical error (right) with mesh size $h=1/64$.

For the CH equation, phase coarsening stage is slow in time, and so one
chooses relatively larger time step to reduce the total amount of calculation.
We choose $\delta=1\times 10^{-3}$ in our experiment that falls under
$\delta>\frac{4\epsilon^{2}}{c^{4}}$. In this case, the iteration count for
the NN method are shown in Table 2.

$h$ $\theta$ |  0.l |  0.2 |  0.25 |  0.3 |  0.4 |  0.5 |  0.6 |  0.7 |  0.8
---|---|---|---|---|---|---|---|---|---
1/64 |  36 |  13 |  2 |  13 |  39 |  - |  - |  - |  -
1/128 |  37 |  13 |  2 |  13 |  40 |  - |  - |  - |  -
1/256 |  38 |  13 |  2 |  13 |  41 |  - |  - |  - |  -
1/512 |  38 |  13 |  2 |  14 |  41 |  - |  - |  - |  -
Table 2: Number of iteration compared for unequal subdomain in 1D.

In Fig 2, on the left we plot the error curves for NN method for different
values of $\theta$ for CH equation in 2D, and on the right we show the
numerical error and theoretical estimates for $\theta=1/4$.

(a)

(b)

Figure 2: Iteration compared in 2D (left), and comparison of theoretical and
numerical error (right) with mesh size $h=1/64$.

##  4 Conclusions

We studied the Neumann-Neumann method for the CH equation for two subdomains.
We proved convergence estimates for the case of one dimensional NN and also
extended the analysis to the 2D CH equation using Fourier techniques. Lastly
we validate our theoretical findings with numerical experiments.

#### Acknowledgement:

I wish to express my appreciation to Dr. Bankim C. Mandal for his constant
support and stimulating suggestions and also like to thank the CSIR India for
the financial assistance and IIT Bhubaneswar for research facility.

## References

  * [1] Bertozzi, A.L., Esedoḡlu, S., Gillette, A.: Inpainting of binary images using the Cahn-Hilliard equation.  IEEE Trans. Image Process. 16(1), 285–291 (2007)
  * [2] Bjørstad, P.E., Widlund, O.B.: Iterative methods for the solution of elliptic problems on regions partitioned into substructures.  SIAM J. Numer. Anal. 23(6), 1097–1120 (1986).  10.1137/0723075
  * [3] Bourgat, J.F., Glowinski, R., Le Tallec, P., Vidrascu, M.: Variational formulation and algorithm for trace operation in domain decomposition calculations.  Ph.D. thesis, INRIA (1988)
  * [4] Cahn, J.W.: On spinodal decomposition.  Acta Metall 9(9), 795–801 (1961)
  * [5] Cahn, J.W., Hilliard, W.: Free energy of a nonuniform system. i. interfacial free energy.  J. Chem. Phys. 28(2), 258–267 (1958)
  * [6] Eyre, D.J.: Unconditionally gradient stable time marching the Cahn-Hilliard equation.  In: Computational and mathematical models of microstructural evolution (San Francisco, CA, 1998), _Mater. Res. Soc. Sympos. Proc._ , vol. 529, pp. 39–46. MRS, Warrendale, PA (1998)
  * [7] Eyre, D.J.: An unconditionally stable one-step scheme for gradient systems.  Unpublished article (1998)
  * [8] Le Tallec, P., De Roeck, Y.H., Vidrascu, M.: Domain decomposition methods for large linearly elliptic three-dimensional problems.  Journal of Computational and Applied Mathematics 34(1), 93–117 (1991)
  * [9] Lee, S., Lee, C., Lee, H.G., Kim, J.: Comparison of different numerical schemes for the cahn-hilliard equation.  J. KSIAM 17(3), 197–207 (2013)
  * [10] Lions, P.L.: On the Schwarz alternating method. I.  In: First International Symposium on Domain Decomposition Methods for Partial Differential Equations (Paris, 1987), pp. 1–42. SIAM, Philadelphia, PA (1988)
  * [11] Toselli, A., Widlund, O.B.: Domain Decomposition Methods, Algorithms and Theory, vol. 34.  Springer-Verlag, Berlin (2005)