#  Dynamical behavior of a nonlocal Fokker-Planck equation for a stochastic
system with tempered stable noise

Li Lin 1, Jinqiao Duan 2, Xiao Wang3 and Yanjie Zhang4

1 Center for Mathematical Sciences, Huazhong University of Science and
Technology
Wuhan 430074, China
email<EMAIL_ADDRESS>
2 Department of Applied Mathematics, Illinois Institute of Technology
Chicago, IL 60616, USA
<EMAIL_ADDRESS>
3 School of Mathematics and Statistics, Henan University
Kaifeng 475001, China
email<EMAIL_ADDRESS>
4 School of Mathematics, South China University of Technology
Guangzhou 510000, China
<EMAIL_ADDRESS>

###### Abstract

We characterize a stochastic dynamical system with tempered stable noise, by
examining its probability density evolution. This probability density function
satisfies a nonlocal Fokker-Planck equation. First, we prove a superposition
principle that the probability measure-valued solution to this nonlocal
Fokker-Planck equation is equivalent to the martingale solution composed with
the inverse stochastic flow. This result together with a Schauder estimate
leads to the existence and uniqueness of strong solution for the nonlocal
Fokker-Planck equation. Second, we devise a convergent finite difference
method to simulate the probability density function by solving the nonlocal
Fokker-Planck equation. Finally, we apply our aforementioned theoretical and
numerical results to a nonlinear filtering system by simulating a nonlocal
Zakai equation.

_Key words:_ Nonlocal Fokker-Planck equation, superposition principle,
Schauder estimate, data assimilation, Zakai equation.

Stochastic dynamic systems are widely used to describe various complex
phenomena. The random fluctuations in complex phenomena are usually non-
Gaussian. How to capture the uncertainly propagates and evolves for the
nonlinear stochastic dynamical systems is an important issue. A popular method
is to obtain the probability density function of the solution paths, which
contains the complete statistical information. In this paper, we prove a
superposition principle that the probability measure-valued solution to this
nonlocal Fokker-Planck equation is equivalent to the martingale solution.
Meanwhile, we devise a convergent finite difference method to simulate the
probability density function by solving the nonlocal Fokker-Planck equation,
which helps us to detect the dynamic behavior of the stochastic dynamic
system. Further, we apply our aforementioned theoretical and numerical results
to a nonlinear filtering system by simulating a nonlocal Zakai equation. The
results established in this paper can be used to examine dynamical behaviors
for financial markets, climate dynamics and physics.

##  1 Introduction

Lévy processes have been observed experimentally in fluid dynamics and
polymers and have been used to describe subrecoil laser cooling, turbulent
fluids, very stiff polymers, and the spectral random walk of a single molecule
embedded in a solid [1, 2, 3, 4, 5]. However, in physical systems, the
variance of any stationary processes is finite. Then an unavoidable cutoff is
always present. For example, in the case of a single molecule embedded in a
solid, due to the minimal length between the molecule and the nearest two-
level systems, a cutoff is present in the distribution of the jumps of the
resonance frequency [6, 7]. Many scholars have introduced a smooth exponential
regression towards zero. This makes it possible to derive an analytic
expression for characteristic function and enables one to replace simulations
by more straightforward calculations. This is the tempered stable processes
what we mention in our paper. Unlike the $\alpha$-stable counterpart, a
tempered stable Lévy process has mean, variance, and moments of all order [8].

The Fokker-Planck equation (FPE) is an important deterministic tool for
quantifying stochastic differential equation. The solution of FPE describes
the evolution of the transition probability density for a stochastic system.
The nonlocal FPE for systems with non-Gaussian noise have been derived
(e.g.,[9, 10]). It is hardly possible to have analytical solution for the
nonlocal FPE except the special drift term [11]. A few authors have considered
numerical simulations of such nonlocal FPE. Cont and Voltchkova [12] presented
a finite difference method for solving parabolic integro-differential
equations with possibly singular kernels, which can be used to price European
and barrier options via stochastic models with tempered Lévy processes. Li and
Deng [13] developed a high order difference scheme for a tempered fractional
diffusion equation on a bounded domain, together with stability and
convergence analysis. Gao et al. [14] developed a fast finite difference
scheme to simulate the nonlocal FPE on either a bounded or infinite domain. Xu
et al. [15, 16] developed the path integral method to solve one-dimensional
space fractional Fokker-Planck-Kolmogorov equations.

Our study is divided into two parts. In the first part, we study the
superposition principle and global parabolic Schauder estimate for the
nonlocal FPE under the natural condition $\alpha+\beta>1$. The operator
$\mathcal{L}$ is not a stable-like operator, i.e., we can not find two
positive constants $c_{1}$ and $c_{2}$ such that

$c_{1}\leq\frac{1}{e^{\lambda|y|}}\leq c_{2},~{}~{}\forall
y\in\mathbb{R}^{d}.$ (1.1)

So we can not apply the recent works of Zhang et al. [17] and Chen et al. [18]
who address the parabolic Dirichlet problem through using a probabilistic
argument. A natural question is whether the global parabolic Schauder estimate
$\|p(x,t)\|_{L^{\infty}\left([0,T],C^{\alpha+\beta}_{b}\right)}\leq
C\|g\|_{C^{1+\beta}_{b}}$ holds. We will answer this in Section $2$. Further,
we construct a finite difference scheme to simulate the nonlocal Fokker-Planck
equation. Under a specified condition, the semi-discrete scheme is shown to
satisfy the discrete maximum principle and to be convergent. In the second
part, we apply the aforementioned results to obtain the strong form of Zakai
equation and conduct its numerical similation.

This paper is organized as follows. In section 2, we study the well-poseness
for the nonlocal FPE associated with stochastic differential equation driven
by tempered stable noise. Further, we construct a convergent finite difference
scheme to simulate the nonlocal FPE, and a numerical experiment is conducted
to confirm the theoretical results. In section 3, we apply the aforementioned
results to nonlinear filtering problem. Some concluding remarks are made in
Section 4.

##  2 METHODOLOGY

A: The nonlocal FPE
In this paper, we consider the following stochastic differential equation
(SDE)

$\mathrm{d}X_{t}=f(X_{t})\mathrm{d}t+\mathrm{d}L_{t},$ (2.2)

where $f$ is a $d$-dimensional Borel measurable function, and $L_{t}$ is a
$d$-dimensional tempered stable Lévy process with triplet $(0,0,\nu)$. This
triplet means that the process has shift zero, diffusion zero, and Lévy jump
measure $\nu$. The jump measure $\nu$ is obtained through multiplying the
$\alpha$-stable Lévy measure by an exponentially decaying function. This
transform is called exponential tilting of the Lévy measure for
$\alpha$-stable Lévy measure (see [19]), i.e.,

$\nu({\rm d}y)=\frac{C_{d,\alpha}{\rm d}y}{e^{\lambda|y|}|y|^{d+\alpha}}.$
(2.3)

Here $C_{d,\alpha}$ is a constant, and $\lambda$ is the positive tempering
parameter (see [20]). The parameter $\alpha\in(0,1)\cup(1,2)$ is called the
stable index.

The corresponding FPE for the SDE (2.2) is the following nonlocal parabolic
equation, i.e.,

$\frac{\partial p}{\partial
t}=-\nabla\cdot\left(f(x)p(x,t)\right)+\mathcal{L}p(x,t),$ (2.4)

where

$\displaystyle\mathcal{L}p(x,t)$
$\displaystyle=\int_{\mathbb{R}^{d}\setminus\\{0\\}}\left[p(x+y,t)-p(x,t)-1_{\\{|y|<1\\}}(y)y\cdot\nabla
p(x,t)\right]\nu(\mathrm{d}y)$ (2.5)
$\displaystyle={\bf{p.v.}}\int_{\mathbb{R}^{d}\setminus\\{0\\}}\left[p(x+y,t)-p(x,t)\right]\nu(\mathrm{d}y).$

In the following, we are interested in studying the dynamical properties and
numerical analysis for the equation (2.4) in the principle value sense.

B: Examing the solution of nonlocal FPE

Let $\mu_{t}(dx)=p(x,t)dx$ be the marginal law of $X_{t}$. By Itô’s formula,
$\mu_{t}$ solves the following Fokker-Planck equation in the distributional
sense, i.e.,

$\partial_{t}\mu_{t}=-\nabla\cdot\left(f(x)\mu_{t}\right)+\int_{\mathbb{R}^{d}\setminus\\{0\\}}\left[\mu_{t}(x+y)-\mu_{t}(x)-1_{\\{|y|<1\\}}(y)y\cdot\frac{\nabla\mu_{t}}{\partial
x}\right]\nu(\mathrm{d}y).$ (2.6)

The proof of the existence and uniqueness of the equation (2.11) is based on
two ingredients. The first ingredient will be to establish the superposition
principle for probability measure-valued solution to a nonlocal FPE. The
second ingredient will then be to establish the global Schauder estimate for
the solution of (2.11).
B1: Superposition Principle
Next we will establish the superposition principle for probability measure-
valued solution to the nonlocal FPE. Throughout this paper, we make the
following assumptions.
Hypothesis H.1. The drift coefficient $f$ is locally bounded, measurable and
there exists a positive constant $C$ such that

$\sup_{x}\langle x,f(x)\rangle^{+}\leq C(1+|x|^{2}),$ (2.7)

where $\langle x,f(x)\rangle^{+}$ denote the non-negative part of $\langle
x,f(x)\rangle.$
Hypothesis H.2. The drift coefficient $f$ is locally $\beta$-Hölder
continuous, $\beta\in(0,1)$, i.e., there exists a positive constant $K_{0}$
such that

$|f(x)-f(y)|\leq
K_{0}|x-y|^{\beta},~{}~{}x,y\in\mathbb{R}^{d},~{}s.t~{}~{}|x-y|\leq 1.$ (2.8)

######  Theorem 2.1 (Superposition principle).

Under Hypothesis $\bf{H.1}$, for every weak solution $(\mu_{t})_{t\geq 0}$ of
nonlocal FPE (2.6), there is a martingale solution
$\mathbb{P}\in\mathcal{M}^{\mu_{0}}_{0}(\widetilde{\mathcal{L}})$ such that

$\mu_{t}=\mathbb{P}\circ X^{-1}_{t},~{}~{}\forall t\geq 0,$ (2.9)

###### Proof.

Using the similar technique as [21, Theorem 1.5], we only need to verify

$\hbar^{\nu}(x):=\int_{|z|>1+|x|}\log\left(1+\frac{|z|}{1+|x|}\right)\nu(dz)<\infty.$

In fact, $\forall\lambda>0$, $\alpha\in(0,1)\cup(1,2)$ and
$\beta\in(0,\alpha\wedge 1)$, there exists a positive constant $C_{\alpha}$,
such that

$\displaystyle\hbar^{\nu}(x)$
$\displaystyle=C_{\alpha}\int_{|z|>1+|x|}\log\left(1+\frac{|z|}{1+|x|}\right)\frac{1}{e^{\lambda|z|}|z|^{1+\alpha}}dz$
(2.10) $\displaystyle\leq
C_{\alpha}\int_{|z|>1+|x|}\log\left(1+\frac{|z|}{1+|x|}\right)\frac{1}{|z|^{1+\alpha}}dz$
$\displaystyle\leq
C_{\alpha}\int_{|z|>1+|x|}\left(\frac{|z|}{1+|x|}\right)^{\beta}\frac{1}{|z|^{1+\alpha}}dz$
$\displaystyle\leq C_{\alpha}\int_{|z|>1}\frac{1}{|z|^{1+\alpha-\beta}}dz$
$\displaystyle\leq C_{\alpha}.$

∎

Next we will show the existence and uniqueness of strong solution for the
nonlocal FPE (2.11).
B2: Schauder Estimate
In this subsection, we are interested in the absolutely continuous solution of
the nonlocal FPE (2.6), which satisfies the time-dependent integro-
differential equation (also called nonlocal FPE), i.e.,

$\frac{\partial p}{\partial
t}=-\nabla\cdot\left(f(x)p(x,t)\right)+\mathcal{L}p(x,t),~{}~{}p(x,0)=g(x).$
(2.11)

Superposition principle tells that the existence and uniqueness of the weak
solution for equation (2.11). To obtain the existence and uniqueness of the
strong solution, we need to establish the global Schauder estimate for the
solution of (2.11).

######  Theorem 2.2.

(Schauder estimate.) Under Hypothesis $\bf{H.1}$ and Hypothesis $\bf{H.2}$.
Suppose that $\operatorname{div}{f}\in
C^{\beta}_{b}(\mathbb{R}^{d},\mathbb{R}),g\in
C^{(1\vee\alpha)+\beta}_{b}(\mathbb{R}^{d},\mathbb{R}),$ and that there exists
a positive constant $C$, such that for the solution $p(x,t)$ of nonlocal FPE
(2.11) with $\alpha\in(0,1)\cup(1,2)$, it holds that

$\|p(x,t)\|_{L^{\infty}\left([0,T],C^{\alpha+\beta}_{b}\right)}\leq
C\|g\|_{C^{(1\vee\alpha)+\beta}_{b}}.$ (2.12)

###### Proof.

See Appendix I. ∎

Then by Theorem 2.2, we obtain the following result.

######  Corollary 2.1.

(Existence and uniqueness of the solution for the nonlocal FPE.) Let
$\alpha\in(0,1)\cup(1,2)$ be fixed, then there exists a unique solution $u\in
C_{b}^{\alpha+\beta}([0,T],\mathbb{R}^{d})$ to the nonlocal FPE (2.11), which
also satisfies the Schauder estimate (2.12).

###### Proof.

Define a family of linear operators by $T_{\theta}=(1-\theta)T_{0}+\theta
T_{1}$, where

$\begin{split}T_{0}p_{\tau}&=p_{\tau}+\mathcal{L}^{*}p_{\tau}-f\nabla
p_{\tau},\\\ T_{1}p_{\tau}&=p_{\tau}+\mathcal{L}^{*}p_{\tau}-f\nabla
p_{\tau}-\operatorname{div}{f}\cdot p_{\tau}.\end{split}$

Note that $T_{0}$ is the case considered in [23, Theorem 3], which maps
$C_{b}^{\alpha+\beta}([0,T],\mathbb{R}^{d})$ onto
$C_{b}^{\beta}([0,T],\mathbb{R}^{d})$. The solvability of the nonlocal FPE
(2.11) is equivalent to the invertibility of the operator $T_{\theta}$. By the
method of continuity, we obtain the required result. ∎

C: Numerical analysis

In this subsection, we devised a convergent finite difference method to
simulate the nonlocal FPE (2.11).
C1: A convergent finite difference scheme
Consider the problem in two cases: one case is on a bounded standard domain
$D=(-1,1)$ with the absorbing condition, and the other is on the infinite
domain $\mathbb{R}$. The absorbing condition is that the “partical” $X_{t}$
disappears or is killed when $X_{t}$ is outside a bounded domain $D$. So the
probability density function $p(x,t)$ of being outside of the bounded domain
$D$ is zero, i.e.,

$p(x,t)=0,\;\;\;x\notin(-1,1).$ (2.13)

In fact, by the following transformation

$\displaystyle\tilde{x}$
$\displaystyle=\frac{2x-(a+b)}{b-a},~{}~{}\tilde{y}=\frac{2y}{b-a},$ (2.14)
$\displaystyle u(\tilde{x},t)$
$\displaystyle=p\left(\frac{(b-a)\tilde{x}+(a+b)}{2},t\right).$

We can also convert the finite interval $\widetilde{{D}}=(a,b)$ into standard
domain ${D}=(-1,1)$. Then the nonlocal FPE (2.11) can be rewritten as

$\displaystyle\frac{\partial u(\tilde{x},t)}{\partial
t}=-\frac{\partial\left(\tilde{f}(\tilde{x})u(\tilde{x},t)\right)}{\partial{{\tilde{x}}}}$
(2.15)
$\displaystyle+\int_{\mathbb{R}\setminus\\{0\\}}\left[u(\tilde{x}+\tilde{y},t)-u(\tilde{x},t)-I_{\\{|\tilde{y}|<\frac{2}{b-a}\\}}(\tilde{y})\tilde{y}\frac{\partial
u(\tilde{x},t)}{\partial\tilde{x}}\right]\frac{\widetilde{C}_{\alpha}}{e^{\tilde{\lambda}|\tilde{y}|}|\tilde{y}|^{1+\alpha}}d\tilde{y},~{}~{}\tilde{x}\in(-1,1),$

where

$\displaystyle\tilde{f}(\tilde{x})$
$\displaystyle=\frac{2}{b-a}f\left(\frac{(b-a)\tilde{x}+a+b}{2}\right),$
(2.16) $\displaystyle\widetilde{C}_{\alpha}$
$\displaystyle=(\frac{b-a}{2})^{-\alpha}C_{\alpha},~{}~{}\tilde{\lambda}=\frac{b-a}{2}\lambda.$

Next, the detailed procedure of the proposed algorithm will be described as
follows:

$(i)$ Step one is to introduce the following function,

$\int_{s}^{\infty}x^{-\varrho}e^{-x}{\rm
d}x=s^{-\frac{\varrho}{2}}e^{-\frac{s}{2}}W_{-\frac{\varrho}{2},\frac{1-\varrho}{2}}(s),\text{for}\;s>0,$
(2.17)

where $W$ is the Whittaker W function.

The principal value of integral
$\int_{\mathbb{R}\setminus\\{0\\}}1_{\\{|y|<1\\}}(y)y\cdot\frac{\partial
p(x,t)}{\partial x}\nu(\mathrm{d}y)$ vanishes. Thus the integral term of
equation (2.11) becomes

$\displaystyle\int_{\mathbb{R}\setminus\\{0\\}}\left[p(x+y,t)-p(x,t)\right]\nu(\mathrm{d}y)$
(2.18)
$\displaystyle=C_{\alpha}\int_{-\infty}^{-1-x}\frac{p(x+y,t)-p(x,t)}{e^{\lambda|y|}|y|^{1+\alpha}}{\rm
d}y+C_{\alpha}\int_{-1-x}^{1-x}\frac{p(x+y,t)-p(x,t)}{e^{\lambda|y|}|y|^{1+\alpha}}{\rm
d}y+C_{\alpha}\int_{1-x}^{\infty}\frac{p(x+y,t)-p(x,t)}{e^{\lambda|y|}|y|^{1+\alpha}}{\rm
d}y$
$\displaystyle=-C_{\alpha}p(x,t)\left[W_{1}(x)+W_{2}(x)\right]+C_{\alpha}\int_{-1-x}^{1-x}\frac{p(x+y,t)-p(x,t)}{e^{\lambda|y|}|y|^{1+\alpha}}{\rm
d}y,$

where

$\displaystyle W_{1}(x)$
$\displaystyle=\lambda^{\frac{\alpha-1}{2}}(1+x)^{-\frac{\alpha+1}{2}}e^{-\frac{\lambda(1+x)}{2}}W_{-\frac{1+\alpha}{2},-\frac{\alpha}{2}}(\lambda(1+x)),$
(2.19) $\displaystyle W_{2}(x)$
$\displaystyle=\lambda^{\frac{\alpha-1}{2}}(1-x)^{-\frac{\alpha+1}{2}}e^{-\frac{\lambda(1-x)}{2}}W_{-\frac{1+\alpha}{2},-\frac{\alpha}{2}}(\lambda(1-x)).$

$(ii)$ Step two is to divide the interval $[-2,2]$ into $4J$ subintervals and
define $x_{j}=jh$ for $-2J\leq j\leq 2J$ integer, where $h=\frac{1}{J}$. Using
central difference scheme for the first and two derivatives and modifying the
“punched-hole” trapezoidal rule in the nonlocal term, we get the
discretization scheme of (2.11), i.e.,

$\displaystyle\frac{d{P_{j}}}{dt}:=$ $\displaystyle
C_{h}\frac{P_{j+1}-2P_{j}+P_{j-1}}{h^{2}}-[(fP)_{x,j}^{+}+(fP)_{x,j}^{-}]$
(2.20) $\displaystyle-
C_{\alpha}\left[W_{1}+W_{2}\right]P_{j}+C_{\alpha}h\sum\limits_{k=-J-j,k\neq
0}^{J-j}\\!{{}^{\prime\prime}}\frac{P_{j+k}-P_{j}}{e^{\lambda|x_{k}|}|x_{k}|^{1+\alpha}},$

where $C_{h}=-C_{\alpha}\zeta(\alpha-1)h^{2-\alpha}$, $(fP)_{x,j}^{+}$ and
$(fP)_{x,j}^{-}$ are defined as the global Lax-Friedrichs flux splitting [24],
i.e., $(fP)^{\pm}=\frac{1}{2}(fP\pm MP)$ with $M=\max{|f(x)|}$. Here
$\sum^{{}^{\prime\prime}}$ denotes the quantities corresponding to the two end
summation indices are multiplied by $1/2$.

It is noted that for solving the initial value problem of (2.11) on the
infinite domain $\mathbb{R}$, the semi-discrete equation becomes

$\displaystyle\frac{d{P_{j}}}{dt}:=$ $\displaystyle
C_{h}\frac{P_{j+1}-2P_{j}+P_{j-1}}{h^{2}}-[(fP)_{x,j}^{+}+(fP)_{x,j}^{-}]+C_{\alpha}h\sum\limits_{k=-J-j,k\neq
0}^{J-j}\\!{{}^{\prime\prime}}\frac{P_{j+k}-P_{j}}{e^{\lambda|x_{k}|}|x_{k}|^{1+\alpha}},$
(2.21)

where $J=\frac{\widetilde{L}}{h}$ and $\widetilde{L}$ is large enough so that
the results are convergent.

($iii$) Step three is to give the condition of numerical scheme satisfying the
discrete maximum principle. For the absorbing boundary condition and forward
Euler scheme for time derivative, the scheme (2.20) satisfies the discrete
maximum principle with $f=0$, if $\bigtriangleup t$ and $h$ satisfy the
following condition,

$\displaystyle\frac{\bigtriangleup
t}{h^{\alpha}}\leq\frac{1}{2C_{\alpha}[1+\frac{1}{\alpha}-\zeta(\alpha-1)]}.$
(2.22)

($iv$) Step four is to illustrate the convergence analysis. The numerical
solution $P^{n}_{j}$ of (2.21) converges to the analytic solution to (2.11)
for $x_{j}$ in $[-\tilde{L}/2,\tilde{L}/2]$ when the refinement path satisfies
(2.22) and the length of the integration interval $2\tilde{L}$ in (2.21) tends
to $\infty$.

######  Remark 2.1.

The detailed proofs of Step ($iii$) and Step ($iv$) see Appendix II and III
respectively.

C2: Numerical experiments
Here we present an example to illustrate our numerical method. We take the
initial condition $p(x,0)=\sqrt{\frac{40}{\pi}}e^{-40x^{2}}$ and the finite
interval $D=(-4,4)$. First, we consider the effect of stability index $\alpha$
in Fig.1. We take $\alpha=0.2,0.6,1.2,1.6$ without drift coefficient $f=0$ at
time $t=0.5$. As we see, the larger the $\alpha$ becomes, the flatter the
probability density function is. Then, we illustrate the evolution for the
probability density function with $\alpha=1.5,f=x-x^{3}$ at different time
$t=0.4,0.8,1.6$ and other factors fixed in Fig.2. As time goes on, the
‘particles’ are gradually tend to stay at the two stable points. The maxima of
the probability density function approach $x=\pm 1$. Next, we present the
effect of tempering parameter $\lambda$ in Fig.3. From the figure, we see that
value of the density become larger near the origin when $\lambda$ is larger.
It is the opposite far away from the origin. We also present the effect of
drift terms ($f=0,x-x^{3}$) for probability density function with
$\alpha=1.5,\lambda=0.01$ at time $t=1$. The particles centered at $x=0$ for
$f=0$. While, for $f=x-x.^{3}$, the points $x=+-1$ are the two stable steady
states. So the particles centered at one of the stable states. Lastly, we add
the Monte Carlo solutions of stochastic differential equations to verify the
correctness of the finite difference method for $\alpha=0.5,\lambda=0.01,f=0$
at time $t=4$ in Fig.5.

Figure 1:  The effect of stability index $\alpha=0.2,0.6,1.2,1.6$ for
probability density function without drift coefficient ($f=0$) and
$\lambda=0.01$ at time $t=0.5$.  Figure 2:  The evolution for probability
density function with drift coefficient $f=x-x^{3},\alpha=1.5,\lambda=0.01$ at
time $t=0.4,0.8,1.6$.  Figure 3:  The effect of positive tempering parameter
$\lambda=0.01,0.05,0.1$ for probability density function without drift
coefficient ($f=0$) at time $t=1$.  Figure 4:  The effect of drift terms
($f=0,x-x^{3}$) for probability density function with
$\alpha=1.5,\lambda=0.01$ at time $t=1$.  Figure 5:  Compare the Monte Carlo
simulation with finite difference method for $\alpha=0.5,\lambda=0.01,f=0$ at
time $t=4$.

##  3  Application to a nonlinear filtering problem

The data assimilation method represented by nonlinear filtering has been
widely used in many fields, which is a procedure to extract system state
information by observations. Here we apply our aforementioned theoretical and
numerical results to a nonlinear filtering system by simulating a nonlocal
Zakai equation.

###  3.1 The strong form of Zakai equation

Consider the following signal-observation systems on $\mathbb{R}^{2}$

$\left\\{\begin{aligned} dX_{t}&=f(X_{t})dt+dL_{t},\\\
dY_{t}&=h(X_{t})dt+d{W_{t}},\end{aligned}\right.$ (3.23)

where $f$ is given deterministic measurable function and $h$ is a bounded
measurable continuous function. $W_{t}$ is the standard Brownian motion, which
is independent of $L_{t}$.

Let

$\mathcal{Z}_{t}=\bf{{\sigma}}(Y_{s}:0\leq s\leq t)\vee\mathcal{N},$ (3.24)

where $\mathcal{N}$ is the collection of all $\mathbb{P}$ -negligible sets of
$(\Omega,\mathcal{F})$. Define
$\mathcal{Z}=\sigma(\bigcup_{t}\mathcal{Z}_{t})$. By the version of Girsanov’s
change of measure theorem, we obtain a new probability measure
$\widetilde{\mathbb{P}}$, such that the observation $Z_{t}$ becomes
$\widetilde{\mathbb{P}}$-independent of the signal variables $(X_{t},Y_{t})$.
This can be done through

$\frac{d\widetilde{\mathbb{P}}}{d\mathbb{P}}=\exp\left(-\sum\limits_{i=1}^{m}{\int_{0}^{t}{{h^{i}}({X_{s}})dW_{s}^{i}}}-\frac{1}{2}\sum\limits_{i=1}^{m}{\int_{0}^{t}{{h^{i}}{{({X_{s}})}^{2}}}}ds\right).$
(3.25)

For every bounded differentiable function $\varphi$, by the Kallianpur-
Striebel formula, we have the following representation

$\mathbb{E}[\varphi(X_{t},Y_{t})|\mathcal{Z}_{t}]=\frac{\widetilde{\mathbb{E}}\left[\widetilde{\mathbb{R}}_{t}\varphi(X_{t},Y_{t})|\mathcal{Z}\right]}{\widetilde{\mathbb{E}}\left[\widetilde{\mathbb{R}}_{t}|\mathcal{Z}\right]},$
(3.26)

where

$\widetilde{\mathbb{R}}_{t}=\left.\frac{d\widetilde{\mathbb{P}}}{d\mathbb{P}}\right|_{\mathcal{Z}_{t}}.$
(3.27)

The unnormalized conditional distribution of $\varphi(X_{t})$, given $Y_{t}$,
is defined as

${{P}_{t}}(\varphi)=\widetilde{\mathbb{E}}\left[\widetilde{\mathbb{R}}_{t}\varphi(X_{t},Y_{t})|\mathcal{Z}_{t}\right].$
(3.28)

Heuristically, if the unconditional distribution of the signal
${P}_{t}(\varphi)$ has a density $p(x,t)$ with respect to Lebesgue measure for
all $t>0$, i.e.,

${P}_{t}(\varphi)=\int_{\mathbb{R}}\varphi(x)p(x,t)dx,$ (3.29)

then the unnormalized density $p(x,t)$ satisfies the following Zakai equation.

######  Theorem 3.1.

Under Hypotheses $\bf{H1}$-$\bf{H2}$ and $\bf{H3}$, the probability density
function $p(x,t)$ satisfies the following Zakai equation, i.e.,

$dp(x,t)={{A}}^{*}p(x,t)dt+h(x)p(x,t)dY_{t},$ (3.30)

where

${{A}}^{*}p(x,t)=-\nabla\cdot\left(f(x)p(x,t)\right)+\int_{{\mathbb{R}^{d}}\backslash\\{0\\}}\left[p(x+z,t)-p(x,t)\right]\nu(dz).$
(3.31)

###### Proof.

The proof is similar to [10, Theorem 5]. ∎

###  3.2  An example

Here we consider the gradient system as the signal system with time-
independent bistable potential, which describes the evolution of physical
phenomena in ideal fluctuating environments.

$\displaystyle dX_{t}$ $\displaystyle=-\frac{\partial V(X_{t},t)}{\partial
x}dt+dL_{t},$ (3.32)

where $V(x,t)=-\frac{1}{2}x^{2}+\frac{1}{4}x^{4}$.

In this experiment, we assume that the observation system is given by

$dY_{t}=\cos(X_{t})dt+dW_{t},$ (3.33)

Using Theorem 3.1, the strong form of Zakai equation for the signal-
observation system (3.32)-(3.33) is

$dp(x,t)={\widetilde{\mathcal{L}}}^{*}p(x,t)dt+\cos x\cdot p(x,t)dY_{t},$
(3.34)

where

$\displaystyle{\widetilde{\mathcal{L}}}^{*}p(x,t)=$ $\displaystyle\
-\frac{\partial\left(f(x)p(x,t)\right)}{\partial
x}+\int_{{\mathbb{R}}\backslash\\{0\\}}\left[p(x+z,t)-p(x,t)\right]\nu(dz).$

In Fig. 6, we simulate the signal-observation processes in Eq. (3.32) and
(3.33)(see [26]). We take $\lambda=0.01,\alpha=1.5,f=x-x^{3},X(0)=-1,Y(0)=-1$
with the terminal time $t=1$. The solution of Zakai equation in Eq. (3.34)
with $p(x,0)=\sqrt{\frac{40}{\pi}}e^{-40x^{2}}$ is illustrated at time $t=1$
in Fig. 7.

Figure 6:  (a) The graph plots the signal process for
$\lambda=0.01,\alpha=1.5,f=x-x^{3},X(0)=-1,Y(0)=-1$. (b) The graph plots the
observation process for $\lambda=0.01,\alpha=1.5,f=x-x^{3},X(0)=-1,Y(0)=-1$.
Figure 7:  A sample path of Zakai equation with
$p(x,0)=\sqrt{\frac{40}{\pi}}e^{-40x^{2}}$ at time $t=1$.

##  4 Conclusion

In this paper, we have studied the superposition principle and numerical
analysis for the nonlocal Fokker-Planck equation associated with a stochastic
dynamical system with tempered stable Lévy noise. Firstly, we have shown the
superposition principle for probability measure-valued solution to a nonlocal
Fokker Planck equation, which yields the equivalence between martingale
problem for the stochastic dynamical system with tempered stable noise and the
corresponding nonlocal parabolic equation. Secondly, we have derived a global
parabolic Schauder estimate for the nonlocal Fokker-Planck equation under the
condition $\alpha+\beta>1$. This leads to the existence and uniqueness have
been obtained in the
$L^{\infty}\left(\left[0,T\right],C^{\alpha+\beta}_{b}\right)$. Thirdly, we
have constructed a finite difference scheme to simulate the nonlocal Fokker-
Planck equation. Under a specific condition, the scheme has been shown to
satisfy the discrete maximum principle and to be convergent. Finally, we have
applied the nonlocal Fokker-Planck equation to obtain the strong form of Zakai
equation and conduct its numerical assimilation. The results established in
this paper can be used to examine dynamical behaviors for financial markets,
climate dynamics and physics.

## Acknowledgements

The research of L.Lin was supported by the NSFC grants 11531006 and 11771449.
The research of X. Wang was supported by the NSFC grant 11901159. The research
of Y. Zhang was supported by the NSFC grant 11901202.

## Appendix

### Appendix I

###### Proof.

Step 1. For the case $\alpha\in(0,1),$ we set $\tau=T-t$. Then the equation
(2.11) can be rewritten as follows

$\begin{cases}p_{\tau}-f\nabla p+\mathcal{L}^{*}p=\operatorname{div}{f}\cdot
p\\\ p(x,T)=g(x),\end{cases}$ (A1)

where $T>0$ is a fixed final time.

In the following, we will decompose the Lévy measure and give a bound for the
character function. Since $\nu(dy)=e^{-\lambda|y|}\cdot\nu_{\alpha}(dy),$
where $\nu_{\alpha}$ is a symmetric stable Lévy measure. By the polar
coordinates $y=r\xi,(r,\xi)\in\mathbb{R_{+}}\times\mathbb{S}^{d-1}$, where
$\mathbb{S}^{d-1}$ represents the $d-1$-dimensional sphere, then the jump
measure $\nu_{\alpha}$ can be decomposed as

$\nu_{\alpha}(B)=\int_{\mathbb{S}^{d-1}}\hat{\mu}(d\xi)\int_{0}^{\infty}1_{B}(r\xi)\frac{dr}{r^{1+\alpha}},\;~{}~{}\text{for}\;B\in\mathcal{B}{(\mathbb{R}^{d})},$

where $\hat{\mu}$ is a finite measure on $\mathbb{S}^{d-1}.$

Note that $\hat{\mu}$ is non-degenerate. Then there exists $\eta\geq 1$, such
that for all $m\in\mathbb{R}^{d},$

$\eta^{-1}|m|^{\alpha}\leq\int_{\mathbb{S}^{d-1}}|(m,\xi)|^{\alpha}\hat{\mu}(d\xi)\leq\eta|m|^{\alpha},~{}~{}\alpha\in(0,1).$
(A2)

The Lévy symbol associated with $\mathcal{L}$ is given

$\begin{split}\psi(\lambda)&=\exp{\left\\{\int_{\mathbb{R}^{d}}\left(e^{(im,y)}-1\right)\cdot
e^{-\lambda|y|}\nu_{\alpha}(dy)\right\\}}\\\
&=\exp{\left\\{\int_{\mathbb{S}^{d-1}}\hat{\mu}(d\xi)\int_{0}^{\infty}(e^{i(m,r\xi)}-1)e^{-\lambda|r\xi|}\frac{dr}{r^{1+\alpha}}\right\\}}\\\
&=\exp{\left\\{\int_{\mathbb{S}^{d-1}}\hat{\mu}(d\xi)\int_{0}^{\infty}\left[\left(e^{(i(m,\xi)-\lambda)r}-1\right)-\left(e^{-\lambda
r}-1\right)\right]\frac{dr}{r^{1+\alpha}}\right\\}}\\\
&:=\exp\\{I_{1}+I_{2}\\}.\end{split}$

Set $\hat{\mu}_{1}=-\Gamma(-\alpha)\hat{\mu}$. Then we have

$\displaystyle I_{1}$
$\displaystyle=\int_{\mathbb{S}^{d-1}}\hat{\mu}(d\xi)\int_{0}^{\infty}\left(e^{(i(m,\xi)-\lambda)r}-1\right)\frac{dr}{r^{1+\alpha}},$
(A3)
$\displaystyle=-\int_{\mathbb{S}^{d-1}}(\lambda-i(m,\xi))^{\alpha}\hat{\mu}_{1}(d\xi)$
(A4)
$\displaystyle=-\int_{\mathbb{S}^{d-1}}|\lambda-i(m,\xi)|^{\alpha}e^{i\alpha\cdot
arg(\lambda-i(m,\xi))}\hat{\mu}_{1}(d\xi),$ (A5)

and

$\displaystyle I_{2}$
$\displaystyle=\int_{\mathbb{S}^{d-1}}\hat{\mu}(d\xi)\int_{0}^{\infty}\left(e^{-\lambda
r}-1\right)\frac{dr}{r^{1+\alpha}}$
$\displaystyle=\lambda^{\alpha}\hat{\mu}_{1}(\mathbb{S}^{d-1}).$ (A6)

Define $\theta=\arg(\lambda-i(m,\xi))$. Then we have

$e^{i\alpha\cdot arg(\lambda-
im\xi)}=e^{i\alpha\theta}=\cos(\alpha\theta)+i\sin(\alpha\theta).$

Since $\hat{\mu}_{1}$ is a symmetric measure, we conclude that

$I_{1}=-\int_{\mathbb{S}^{d-1}}\left((m,\xi)^{2}+\lambda^{2}\right)^{\alpha/2}cos(\alpha\theta)\hat{\mu}_{1}(d\xi).$

Recall that $I_{1}$ is the Lévy symbol, $\mathscr{R}I_{1}\leq 0$. There exists
$\eta>1,$ such that

$-\eta|m|^{\alpha}-C\leq I_{1}\leq-\eta^{-1}|m|^{\alpha},$ (A7)

where $C$ is a positive constant.

Combined the inequality (A7) with equality (A6), we have

$-\eta|m|^{\alpha}+\lambda^{\alpha}\hat{\mu}_{1}(\mathbb{S}^{d-1})-C\leq
I_{1}+I_{2}\leq\lambda^{\alpha}\hat{\mu}_{1}(\mathbb{S}^{d-1})-\eta^{-1}|m|^{\alpha}$
(A8)

Step 2. We will establish the upper bound for the heat-kernel and for the
derivatives of the heat-kernel. In fact, for the operator $\mathcal{L}$, it is
well-known that the associated convolution Markov semigroup has a $C^{\infty}$
density $p_{\alpha}(\cdot,t)$. By Fourier inversion, for all $t>0,$
$y\in\mathbb{R}^{d}$, we get

$p_{\alpha}(y,t)=\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\exp(-i(y,m))\exp\left(t(I_{1}+I_{2}\right)dm.$

Moreover, we have the following two results :

  1. 1.

For some $c=c(\gamma,\alpha)>0$ and $t\in(0,1]$, there exists $\alpha\in(0,1)$
such that for every $\gamma$ satisfying $0<\gamma<\alpha$, we have
$\int_{\mathbb{R}^{d}}|y|^{\gamma}p_{\alpha}(y,t)dy\leq ct^{\gamma/{\alpha}}.$

  2. 2.

There exists a positive constant $c=c(\alpha,\beta)>0$ such that
$\int_{\mathbb{R}^{d}}|y|^{\beta}|D_{y}^{k}p_{\alpha}(y,t)|dy\leq\frac{c}{t^{[k-\beta]/{\alpha}}},\;t\in(0,1],k=1,2.$
Here $D_{y}^{1}p_{\alpha}(y,t)$ and $D_{y}^{2}p_{\alpha}(y,t)$ denote the
gradient and Hessian matrix in the $y$-variable, respectively.

On the one hand, by the inequality (A8), we have

$e^{(I_{1}+I_{2})}\leq
e^{\lambda^{\alpha}\hat{\mu}_{1}(\mathbb{S}^{d-1})}e^{-\eta^{-1}|m|^{\alpha}}\leq
Ce^{-\eta^{-1}|m|^{\alpha}}.$

Thus

$\begin{split}\int_{\mathbb{R}^{d}}|y|^{\gamma}p_{\alpha}(y,t)dy&\leq
C\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}e^{-i(y,m)}|y|^{\gamma}e^{-t\eta^{-1}|m|^{\alpha}}dmdy\\\
&\leq
C\int_{\mathbb{R}^{d}}\left(\int_{\mathbb{R}^{d}}e^{-i(y,m)}|y|^{\gamma}dy\right)e^{-t\eta^{-1}|m|^{\alpha}}dm\\\
&\leq C\int_{\mathbb{R}^{d}}|m|^{-\gamma-1}e^{-t\eta^{-1}|m|^{\alpha}}dm\\\
&\leq
C\int_{\mathbb{R}^{d}}e^{-t\eta^{-1}|m|^{\alpha}}(t\eta^{-1}|m|^{\alpha})^{\frac{-\gamma-\alpha}{\alpha}}t^{\frac{\gamma+\alpha}{\alpha}}\frac{1}{t}d(-t\eta^{-1}|m|^{\alpha})\\\
&\leq Ct^{\frac{\gamma}{\alpha}}.\end{split}$

On the other hand,

$\begin{split}\int_{\mathbb{R}^{d}}|y|^{\beta}|D_{y}^{k}p_{\alpha}(y,t)|dy&\leq\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}|m|^{k}e^{-i(y,m)}|y|^{\beta}e^{-t\eta^{-1}|m|^{\alpha}}dmdy\\\
&\leq C\int_{\mathbb{R}^{d}}|m|^{k-\gamma-1}e^{-t\eta^{-1}|m|^{\alpha}}dm\\\
&\leq
Ct^{\frac{\beta-k}{\alpha}}\leq\frac{C}{t^{[k-\beta]/{\alpha}}}.\end{split}$

Step 3. Since $\operatorname{div}{f}\in C^{\beta}_{b}$, which will not affect
the discuss as in [23, Theorem 3]. We also have

$\|p(x,t)\|_{L^{\infty}\left([0,T],C^{\alpha+\beta}_{b}\right)}\leq
C\|g\|_{C^{1+\beta}_{b}}.$ (A9)

Step 4. For the case $\alpha\in(1,2),$ the condition $\alpha+\beta>1$ is
obviously satisfied. Further, the characteristic function is the same as the
case of $\alpha\in(0,1).$ We can obtain the upper bound for the heat-kernel
and the derivatives of the heat-kernel in the same way as $\alpha\in(0,1)$,
i.e.,

$\|p(x,t)\|_{L^{\infty}\left([0,T],C^{\alpha+\beta}_{b}\right)}\leq
C\|g\|_{C^{\alpha+\beta}_{b}}.$ (A10)

Therefore, for the case $\alpha\in(0,1)\cup(1,2)$, we have

$\|p(x,t)\|_{L^{\infty}\left([0,T],C^{\alpha+\beta}_{b}\right)}\leq
C\|g\|_{C^{(1\vee\alpha)+\beta}_{b}}.$ (A11)

Moreover, the global Schauder estimates still holds. ∎

### Appendix II

###### Proof.

Let $M$ be the maximum value of the initial probability density, then for the
numerical solution $0<P_{j}^{n}\leq M$, we have $\zeta(\alpha-1)\leq 0$.
Applying the explicit Euler to the semi-discrete scheme (2.20), we get

$\displaystyle P_{j}^{n+1}$ $\displaystyle=$ $\displaystyle
P_{j}^{n}-\bigtriangleup
tC_{\alpha}\zeta(\alpha-1)h^{-\alpha}(P_{j+1}^{n}-2P_{j}^{n}+P_{j-1}^{n})$
$\displaystyle-C_{\alpha}\bigtriangleup
t\left[W_{1}(x_{j})+W_{2}(x_{j})\right]P_{j}^{n}+C_{\alpha}\bigtriangleup
th\sum\limits_{k=-J-j,k\neq
0}^{J-j}\\!{{}^{\prime\prime}}\frac{P_{j+k}^{n}-P_{j}^{n}}{e^{\lambda|x_{k}|}|x_{k}|^{1+\alpha}}$
$\displaystyle=$ $\displaystyle\bigg{[}1+2\bigtriangleup
tC_{\alpha}\zeta(\alpha-1)h^{-\alpha}-C_{\alpha}(W_{1}(x_{j})+W_{2}(x_{j}))\bigtriangleup
t$ $\displaystyle-C_{\alpha}\bigtriangleup th\sum\limits_{k=-J-j,k\neq
0}^{J-j}\\!{{}^{\prime\prime}}\frac{1}{e^{\lambda|x_{k}|}|x_{k}|^{1+\alpha}}\bigg{]}P_{j}^{n}-\bigtriangleup
tC_{\alpha}\zeta(\alpha-1)h^{-\alpha}(P_{j-1}^{n}+P_{j+1}^{n})$
$\displaystyle+C_{\alpha}\bigtriangleup th\sum\limits_{k=-J-j,k\neq
0}^{J-j}\\!{{}^{\prime\prime}}\frac{P_{j+k}^{n}}{e^{\lambda|x_{k}|}|x_{k}|^{1+\alpha}}.$

Set

$\displaystyle L:$
$\displaystyle=C_{\alpha}(W_{1}(x_{j})+W_{2}(x_{j}))\bigtriangleup
t+C_{\alpha}\bigtriangleup th\sum\limits_{k=-J-j,k\neq
0}^{J-j}\\!{{}^{\prime\prime}}\frac{1}{e^{\lambda|x_{k}|}|x_{k}|^{1+\alpha}}.$
(B1)

Then by the inequality (2.22), we have

$\displaystyle L$ $\displaystyle\leq C_{\alpha}\bigtriangleup
t\left[\int_{-\infty}^{-1-x_{j}}\frac{{\rm d}y}{e^{\lambda
y}|y|^{1+\alpha}}+\int_{1-x_{j}}^{\infty}\frac{{\rm d}y}{e^{\lambda
y}|y|^{1+\alpha}}+\frac{2h}{e^{\lambda
h}h^{1+\alpha}}+\int_{(-1-x_{j}+\frac{h}{2},1-x_{j}-\frac{h}{2})\backslash(-h,h)}\frac{{\rm
d}y}{e^{\lambda y}|y|^{1+\alpha}}\right]$ (B2) $\displaystyle\leq
C_{\alpha}\bigtriangleup t\left[\frac{2}{e^{\lambda
h}h^{\alpha}}+2\int_{h}^{\infty}\frac{{\rm d}y}{e^{\lambda
y}|y|^{1+\alpha}}\right]$ $\displaystyle\leq\frac{2C_{\alpha}\bigtriangleup
t}{h^{\alpha}}(1+\frac{1}{\alpha})$ $\displaystyle\leq 1+2\bigtriangleup
tC_{\alpha}\zeta(\alpha-1)h^{-\alpha},$

where $x_{j}=jh$.

Therefore we have

$P_{j}^{n+1}\leq\left\\{1-C_{\alpha}\bigtriangleup
t[W_{1}(x_{j})+W_{2}(x_{j})]\right\\}M\leq M.$ (B3)

∎

### Appendix III

###### Proof.

Set $e^{n}_{j}=P^{n}_{j}-p(x_{j},t_{n})$, where $p(x_{j},t_{n})$ is the
analytic solution at the point $(x_{j},t_{n})$, then we have

$\displaystyle e^{n+1}_{j}-e^{n}_{j}$ $\displaystyle=-\Delta
tC_{\alpha}\zeta(\alpha-1)h^{2-\alpha}\frac{e^{n}_{j-1}-2e^{n}_{j}+e^{n}_{j+1}}{h^{2}}$
(C1) $\displaystyle+\Delta tC_{\alpha}h\sum^{J-j}_{k=-J-j,k\neq
0}\frac{e_{j+k}^{n}-e_{j}^{n}}{e^{\lambda|x_{k}|}|x_{k}|^{1+\alpha}}-\Delta
tT^{n}_{j},$

where

$\displaystyle
T^{n}_{j}=C_{\alpha}\zeta(\alpha-1)h^{2-\alpha}\frac{p(x_{j-1},t_{n})-2p(x_{j},t_{n})+p(x_{j+1},t_{n})}{h^{2}}$
(C2) $\displaystyle-C_{\alpha}h\sum^{J-j}_{k=-J-j,k\neq
0}\frac{p(x_{j+k},t_{n})-p(x_{j},t_{n})}{e^{\lambda|x_{k}|}|x_{k}|^{1+\alpha}}+\frac{p(x_{j},t_{n+1})-p(x_{j},t_{n})}{\Delta
t}.$

Further, we have

$\displaystyle
T^{n}_{j}=\frac{1}{2}\frac{\partial^{2}p(x_{j},t_{n})}{\partial^{2}t}\Delta
t+A_{\alpha}\zeta(\alpha-1)h^{4-\alpha}\frac{\partial^{4}p(x_{j},t_{n})}{\partial^{4}x}-B_{\alpha}\zeta(\alpha-3)h^{4-\alpha}\frac{\partial^{4}p(x_{j},t_{n})}{\partial^{4}x}$
(C3) $\displaystyle-D_{\alpha}\frac{\partial}{\partial
y}\left(\frac{p(x_{j}+y,t_{n})-p(x_{j},t_{n})}{e^{\lambda|y|}|y|^{1+\alpha}}\right)\Big{|}^{y=L-x_{j}}_{y=-L-x_{j}}$
$\displaystyle+E_{\alpha}\int_{\\{-\infty,-L-x_{j}\\}\bigcup\\{L-x_{j},\infty\\}}\frac{p(x_{j}+y,t_{n})-p(x_{j},t_{n})}{e^{\lambda|y|}|y|^{1+\alpha}}dy+\cdots,$

where $\zeta(\tau)$ is the Riemann zeta function initially defined for
$\mathbb{R}e\tau>1$ by $\zeta(\tau)=\sum_{k=1}^{\infty}k^{-\tau}$,
$A_{\alpha},B_{\alpha},D_{\alpha},E_{\alpha}$ are constants depending on
$\alpha$.

Obviously, we have

$|T^{n}_{j}|\leq O(\Delta t)+O(h^{2})+O(L^{-\alpha}):=\widetilde{{T}}.$ (C4)

Therefore, the truncation error is uniformly bounded. By the condition (2.22),
we have

$\max|e^{n+1}_{j}|\leq\max|e^{n}_{j}|+\Delta t\widetilde{T}\leq n\Delta
t\widetilde{{T}}.$ (C5)

∎

## Data Availability Statements

The data that support the findings of this study are openly available in
GitHub, Ref. [27].

## References

  * [1] B. Lisowski, D. Valenti, B. Spagnolo, M. Bier, E. Gudowska-Nowak, “Stepping molecular motor amid Lévy white noise,” Phys. Rev. E, 91: 042713(2015).
  * [2] T. Srokowski, “Fluctuations in multiplicative systems with jumps,” Phys. Rev. E, 87: 032104(2013).
  * [3] P. D. Ditlevsen, “Observation of $\alpha$-stable noise induced millennial climate changes from an ice record,” Geophys. Res. Lett., 26: 1441-1444(1999).
  * [4] A. V. Milovanov, J. J. Rasmussen, “Lévy flights on a comb and the plasma staircase,” Phys. Rev. E, 98: 022208(2018).
  * [5] G. Samorodnitsky, M. S. Taqqu, Stable non-Gaussian random process (Chapman & Hall/CRC, New York, 1994).
  * [6] J. Rosinski, “Tempering stable processes,” Stoch. Proc. Appl, 117: 677-707(2007).
  * [7] U. Kuchler, S. Tappe, “Exponential stock models driven by tempered stable processes,” J. Econometrics, 181: 53-63(2014).
  * [8] R. Cont, P. Tankov, Financial modelling with Jump processes (Chapman & Hall/CRC, New York, 2004).
  * [9] D. Schertzer, M. Larchevêque, J. Duan, V. V. Yanovsky, S. Lovejoy. “Fractional Fokker-Planck equation for nonlinear stochastic differential equations driven by non-Gaussian Lévy stable noises,” Journal of Mathematical Physics, 42: 200-212(2001).
  * [10] Y. Zhang, X. Wang, Q. Huang, J. Duan, T. Li, “Numerical analysis and applications of Fokker-Planck equations for stochastic dynamical systems with multiplicative $\alpha$-stable noises,” Appl. Math. Model., 87: 711-730(2020).
  * [11] W. Zan, Y. Xu, J. Kurths, A.V. Chechkin, R. Metzler, “Stochastic dynamics driven by combined Lévy-Gaussian noise: fractional Fokker-Planck-Kolmogorov equation and solution,” J. Phys. A-Math. Theor., 53(2020).
  * [12] R. Cont, E. Voltchkova, “A finite difference scheme for option prices in jump diffusion and exponential Lévy models,” SIAM J. Numer. Anal., 43: 1596-1626(2006).
  * [13] C. Li, W. Deng, “High order schemes for the tempered Fractional diffusion equations,” Adv. Comput. Math, 42: 543-572(2016).
  * [14] T. Gao, J. Duan, X. Li, “Fokker-Planck equations for stochastic dynamical systems with symmetric Lévy process,” Appl. Math. Comput., 278: 1-20(2016).
  * [15] Y. Xu, W. Zan, W. Jia, J. Kurths, “Path integral solutions of the governing equation of SDEs excited by Lévy white noise,” J. Comput. Phys., 394: 41-55(2019).
  * [16] W. Zan, Y. Xu, R. Metzler, J. Kurths, “ First-passage problem for stochastic differential equations with combined parametric Gaussian and Lévy white noises via path integral method,” J. Comput. Phys., 435: 110264(2021).
  * [17] X. Zhang, G. Zhao, “Dirichlet problem for supercritical non-local operators,” arXiv: 1809.05712(2018).
  * [18] Z. Chen, X. Zhang, “Heat kernels and analyticity of non-symmetric jump diffusion semigroups,” Probab. Theory Relat. Fields, 165: 267-312(2016).
  * [19] W. Deng, B. Li, W. Tian, P. Zhang, “Boundary problems for the fractional and tempered fractional operators,” Multiscale Model. Simul., 16:125-149(2018).
  * [20] Y. Yan, W. Deng, D. Nie, “A finite difference approximation for the one and two dimensional tempered fractional Laplacian,” Comm. App. Math. Com., 2: 129-145(2020).
  * [21] M. Röckner, L. Xie, X. Zhang, “Superposition principle for nonlocal Fokker-Planck-Kolmogorov operators,” Probab. Theory Relat. Fields, Doi 10.1007/s00440-020-00985-8.
  * [22] D. W. Strook, S. R. S. Varadhan, Multidimensional diffusion processes(Springer, Berlin, 2006).
  * [23] D. Raynal, P. Chaudru , S. Menozzi, E. Priola, “Schauder estimates for drifted fractional operators in the supercritical case,” J. Funct. Anal., (2019).
  * [24] N. Črnjarić-Žic, S. Vuković, L. Sopta, “On different flux splittings and flux functions in WENO schemes for balance laws,” Comput. Fluids, 35: 1074-1092(2006).
  * [25] C. Shu, S. Osher, “Efficient implementation of essentially non-oscillatory shock-capturing schemes,” J. Comput. Phys, 77(2): 439-471(1988).
  * [26] E. Jum, “Numerical approximation of stochastic differential equations driven by Lévy motion with Infinitely Many Jumps,” Doctoral Dissertations(2015).
  * [27] L. Lin, J. Duan, X. Wang, Y. Zhang, “Dynamical behavior of a nonlocal Fokker-Planck equation for a stochastic system with tempered stable noise,” GitHub, https://github.com/Lin518-hust/Lin518-hust..


# Sharp Morrey regularity theory for a fourth order geometrical equation

Chang-Lin Xiang and Gao-Feng Zheng  Three Gorges Mathematical Research Center,
China Three Gorges University, 443002, Yichang, People’s Republic of China
<EMAIL_ADDRESS>School of Mathematics and Statistics, Central China
Normal University, Wuhan 430079, P. R. China<EMAIL_ADDRESS>
###### Abstract.

This paper is a continuation of the recent work of Guo-Xiang-Zheng [10]. We
deduce sharp Morrey regularity theory for weak solutions to the fourth order
nonhomogeneous Lamm-Rivière equation

$\Delta^{2}u=\Delta(V\nabla u)+div(w\nabla u)+(\nabla\omega+F)\cdot\nabla
u+f\qquad\text{in }B^{4},$

under smallest regularity assumptions of $V,w,\omega,F$ and that $f$ belongs
to some Morrey spaces, which was motivated by many geometrical problems such
as the flow of biharmonic mappings. Our results deepens the $L^{p}$ type
regularity theory of [10], and generalizes the work of Du, Kang and Wang [4]
on a second order problem to our fourth order problems.

###### Key words and phrases:

Keywords: fourth order elliptic equation, regularity theory, Morrey space,
Decay estimates, Riesz potential

###### 1991 Mathematics Subject Classification:

2020 Mathematics Subject Classification: 35J47, 35B65

*: corresponding author

The authors are financially supported by the National Natural Science
Foundation of China (No. 12271296, 12271195).

###### Contents

  1. 1 Introduction and main results
  2. 2 Proof of theorem 1.2
  3. 3 Proof of theorem 1.3
  4. A Morrey spaces and Riesz potential

##  1\. Introduction and main results

In this note, we consider the fourth order elliptic problem

(1.1) $\Delta^{2}u=\Delta(V\nabla u)+div(w\nabla
u)+(\nabla\omega+F)\cdot\nabla u+f\qquad\text{in }B^{4},$

where $B^{4}=B_{1}^{4}(0)\subset\mathbb{R}^{4}$ is the unit ball, $u\in
W^{2,2}(B^{4},\mathbb{R}^{m})$ is a weak solution, and $V,w,\omega,F$ satisfy
the smallest regularity assumption

$\displaystyle V\in
W^{1,2}(B^{4},M_{m}\otimes\Lambda^{1}\mathbb{R}^{4}),\quad\quad w\in
L^{2}(B^{4},M_{m}),$ $\displaystyle\omega\in L^{2}(B^{4},so_{m}),\quad\quad
F\in L^{\frac{4}{3},1}(B^{4},M_{m}\otimes\Lambda^{1}\mathbb{R}^{4}),$

which makes the equation critical when discussing regularity theory.

When $f\equiv 0$, the corresponding homogeneous problem

(1.2) $\Delta^{2}u=\Delta(V\nabla u)+div(w\nabla
u)+(\nabla\omega+F)\cdot\nabla u\qquad\text{in }B^{4}$

was first formulated in the interesting work [13] of Lamm and Rivière, which
was intended to extend the celebrated work of Rivière [15] on the second order
problem

(1.3) $-\Delta u=\Omega\cdot\nabla u\qquad\text{in }B^{2};$

and also to provide a new approach to some fourth order conformally invariant
variational problems such as biharmonic mappings (see e.g. [3, 20, 21, 22] for
more regularity theory on biharmonic mappings). The reader is also referred to
[13] for details on how to write the equation of biharmonic mappings into the
form (1.2).

The flow of extrinsic or intrinsic biharmonic mappings also attracted
extensive study, see e.g. Gastel [5], Lamm [12], Wang [19]. In this case, one
is naturally motivated to consider the nonhomogeneous problem (1.1). For
instance, using the advantage of equation (1.1), the existence of global flow
of extrinsic biharmonic mappings was obtained in [13] for any initial data. We
also mention that another slightly different fourth order elliptic system but
with the same nature as (1.2) was formulated by Struwe [18] and further
studied by Guo, Wang and Xiang [7] and many others. The main results of Lamm
and Rivière [13] can be summarized as follows.

######  Theorem 1.1.

([13]) (1) For any $m\in\mathbb{N}$, there exist constants $C_{m}>0$ and
$\epsilon_{m}>0$ such that if

(1.4)
$\|V\|_{W^{1,2}(B^{4})}+\|w\|_{L^{2}(B^{4})}+\|\omega\|_{L^{2}(B^{4})}+\|F\|_{L^{4/3,1}(B^{4})}<\epsilon_{m},$

then there exist $A\in W^{2,2}\cap L^{\infty}(B_{1/2}^{4},M(m))$ and $B\in
W^{1,4/3}(B_{1/2}^{4},M(m)\otimes\wedge^{2}\mathbb{R}^{4})$ with

(1.5) $\|A\|_{W^{2,2}(B_{1/2}^{4})}+\|{\rm
dist}(A,SO_{m})\|_{L^{\infty}(B_{1/2}^{4})}+\|B\|_{W^{1,4/3}(B_{1/2}^{4})}\leq
C_{m}\epsilon_{m}$

such that

$\nabla\Delta A+\Delta AV-\nabla
Aw+A(\nabla\omega+F)=\sum_{l,k=1}^{4}\partial_{l}B_{lk}\partial_{k}.$

(2) Consequently, under the smallness assumption (1.4), $u$ solves problem
(1.2) in $B_{1/2}^{4}$ if and only if $u$ satisfies the conservation law

(1.6) $\Delta(A\Delta u)=\text{\rm div}\,(K)\qquad\text{in }B_{1/2}^{4},$

where

(1.7) $K\equiv 2\nabla A\Delta u-\Delta A\nabla u+Aw\nabla u-\nabla
AV\cdot\nabla u+A\nabla(V\cdot\nabla u)+B\cdot\nabla u.$

(3) Every weak solution of equation (1.2) is continuous.

For applications of this result to the flow of existence biharmonic mappings,
the reader is referred to the same work [13].

As an improvement, Guo and Xiang [8, 9] showed that every weak solution of
equation (1.2) is in fact locally Hölder continuous. Quite recently, Guo,
Xiang and Zheng [10] considered $L^{p}$ type regularity theory of problem
(1.1). They proved that if $f\in L^{p}(B^{4})$ for some $1<p<4/3$, then every
weak solution $u\in W_{\text{\rm loc}}^{3,p_{1}}(B^{4})$, where
$p_{1}=4p/(4-p)$; while in the case $p=1$ and $f\in L\log L(B^{4})$ there
holds $u\in W_{\text{\rm loc}}^{3,4/3,1}(B^{4})$. The result is optimal in the
sense that, the regularity of $u$ can not be improved too much even if $f\in
L^{p}(B^{4})$ for some $p\geq 4/3$. Also, the third order derivative of $u$ is
the best possible: counterexample shows that in general weak solutions may not
have fourth order weak derivatives.

As a counterpart of the second order problem (1.3), we mention that Rivière
[15] was the first one to formulate problem (1.3) in order to study second
order conformally invariant variational problems, and he also established the
method of conservation law and continuity of weak solutions to problem (1.3).
Sharp and Topping [17] considered the nonhomogeneous problem

(1.8) $-\Delta u=\Omega\cdot\nabla u+f\qquad\text{in }B^{2}.$

They proved, among other results, that if $f\in L^{p}(B^{2})$ for some
$1<p<2$, then every weak solution $u\in W^{1,2}(B^{2})$ of problem (1.8)
belongs to $W_{\text{\rm loc}}^{2,p}(B^{2})$; and in the case $p=1$ and $f\in
L\log L(B^{2})$, then $u\in W_{\text{\rm loc}}^{2,1}(B^{2})$. Inspired by the
work of Wang [20] on harmonic mappings, recently Du, Kang and Wang [4] extend
the theory of Sharp and Topping [17] to the case that $f$ belongs to some
Morrey spaces (for definitions of Morrey spaces, see Appendix A). They proved
that if $f\in M^{1,\lambda}(B^{2})$ for some $0<\lambda<1$, then $\nabla u\in
M_{\text{\rm loc}}^{2,2\lambda}(B^{2})$; and if $f\in M^{p,p\lambda}(B^{2})$
for some $1<p<2$ and $0\leq\lambda<(2-p)/p$, then $\nabla^{2}u\in M_{\text{\rm
loc}}^{p,p\lambda}(B^{2})$ and $\nabla u\in M_{\text{\rm
loc}}^{p^{*},p^{*}\lambda}(B^{2})$. This also extends the Hölder regularity
theory of Wang [20] on harmonic mappings to the general second order problem
(1.8), which maybe have potential applications to nonhomogeneous prescribed
mean curvature equations.

This note can be viewed as a continuation of [10]. Our motivation is to
further extend the $L^{p}$ type regularity theory of [10] to the case that $f$
belongs to some Morrey spaces, and thus extend the Morrey regularity theory of
Du, Kang and Wang [4] to the fourth order problem (1.1). Out first result
reads as follows.

######  Theorem 1.2.

Suppose $u\in W^{2,2}(B^{4})$ is a weak solution to equation (1.1) with $f\in
M^{1,\lambda}(B^{4})$ for some $0<\lambda<1$. Then $\nabla u\in M_{\text{\rm
loc}}^{4,4\lambda}(B^{4})$, $\nabla^{2}u\in M_{\text{\rm
loc}}^{2,2\lambda}(B^{4})$ with

(1.9) $\|\nabla
u\|_{M^{4,4\lambda}(B_{1/2}^{4})}+\|\nabla^{2}u\|_{M^{2,2\lambda}(B_{1/2}^{4})}\leq
C\left(\|u\|_{W^{2,2}(B^{4})}+\|f\|_{M^{1,\lambda}(B^{4})}\right).$

As a result, we have $u\in C_{\text{\rm loc}}^{0,\lambda}(B^{4})$ with

(1.10) $\|u\|_{C^{0,\lambda}(B_{1/2}^{4})}\leq
C\left(\|u\|_{W^{2,2}(B^{4})}+\|f\|_{M^{1,\lambda}(B^{4})}\right),$

where $C>0$ is a positive constant depending only on $m,\lambda$ and the
coefficient functions $V,w,\omega,F$.

For the definition of Morrey spaces, see Appendix A. The result is optimal in
the following sense. Consider the simplest case $\Delta^{2}u=f\in
L^{p}(B^{4})$ for some $1<p<4/3$. Then $f\in L^{p}(B^{4})\subset
M^{1,\lambda}(B^{4})$ with $\lambda=4(1-1/p)\in(0,1)$. By the standard
elliptic regularity theory, $u\in W_{\text{\rm loc}}^{4,p}(B^{4})\subset
C_{\text{\rm loc}}^{0,\lambda}(B^{4})$. This shows that the Hölder regularity
(1.10) is optimal. So it is natural to ask whether this result hold for
$\lambda=1$. The answer is negative, even under stronger assumptions on $f$.
Indeed, consider the case $p=4/3$. Then $f\in L^{4/3}(B^{4})\subset
M^{1,1}(B^{4})$. It is known from elliptic regularity theory that the best
possible regularity is $u\in W_{\text{\rm loc}}^{4,4/3}(B^{4})$ and it is
possible that $u\not\in C_{\text{\rm loc}}^{0,1}(B^{4})$.

In the case $f$ has better Morrey regularity, we have the following improved
regularity.

######  Theorem 1.3.

Suppose $u\in W^{2,2}(B^{4})$ is a weak solution to equation (1.1) with $f\in
M^{p,p\lambda}(B_{1}^{4})$ for some $1<p<4/3$ and $0\leq\lambda<(4-3p)/p$.
Then $\nabla u\in M_{\text{\rm loc}}^{p_{3},p_{3}\lambda}(B^{4})$,
$\nabla^{2}u\in M_{\text{\rm loc}}^{p_{2},p_{2}\lambda}(B^{4})$,
$\nabla^{3}u\in M_{\text{\rm loc}}^{p_{1},p_{1}\lambda}(B^{4})$ with
$p_{i}=4p/(4-ip)$ ($1\leq i\leq 3$), and

$\sum_{i=1}^{3}\left\|\nabla^{i}u\right\|_{M^{p_{4-i},p_{4-i}\lambda}(B_{1/2}^{4})}\leq
C\left(\|u\|_{W^{2,2}(B^{4})}+\|f\|_{M^{1,\lambda}(B^{4})}\right)$

holds for some constant $C>0$ depending only on $m,p,\lambda$ and the
coefficient functions $V,w,\omega,F$.

As $M^{p,0}(B^{4})=L^{p}(B^{4})$, this result can be viewed as an extension of
the $L^{p}$ type regularity theory of Guo, Xiang and Zheng [10]. As already
explained in [10], in general, the third order Sobolev regularity is the best
possible. Thus in our case we still have as most third order regularities.
Moreover, the indexes are all optimal, which can be seen from the simple case
$\Delta^{2}u=f$. We leave the details for interested readers.

The strategy of proving Theorems 1.2 and 1.3 can be briefly illustrated as
follows. Theorems 1.2 follows from a Morrey type decay estimates, in which the
conservation law of Theorem 1.1 is crucial. To prove Theorem 1.3, we will
begin with the $W^{3,p_{1}}$-regularity theory of (1.1) in view of [10]. So
under the assumption $f\in M^{p,p\lambda}(B^{4})$, we have $u\in
W^{3,p_{1}}(B_{3/4}^{4})$. Next we use a perturbation method and duality
method to derive the Morrey regularity of $\nabla^{i}u$ ($i=1,2,3$).

Our notations are standard. Throughout we use $A\lesssim B$ to mean that
$A\leq CB$ for some constant $C>0$ depending only $m,p,\lambda$ and
independent of the parameter $\epsilon$. We use $I_{\alpha}=|x|^{\alpha-n}$ to
denote Riesz potential for all $0<\alpha<n$, such that by $I_{\alpha}(f)$ we
mean

$I_{\alpha}(f)=|x|^{\alpha-n}\ast f.$

For Riesz potential theory on Morrey spaces, see Proposition (A.1) in the
Appendix.

##  2\. Proof of theorem 1.2

In this section we prove theorem 1.2 using a Morrey type decay estimates.

Using a standard translation and scaling argument (see [10, subsection 2.3]
for details), we can assume

$\|V\|_{W^{1,2}(B^{4})}+\|w\|_{L^{2}(B^{4})}+\|\omega\|_{L^{2}(B^{4})}+\|F\|_{L^{4/3,1}(B^{4})}<\epsilon,$

where $\epsilon<\epsilon_{m}$ so that by the conservation law of Theorem 1.1
we have

$\Delta{\rm div}(A\nabla u)={\rm div}\tilde{K}+Af\qquad\text{in }B_{1/2}^{4},$

where

$\tilde{K}=K+\Delta(\nabla A\cdot\nabla u)$

with $A,B$ and $K$ being defined as in Theorem 1.1.

To study $A\nabla u$, we use the Hodge decomposition to deduce

$Adu=d\tilde{u}_{1}+d^{\ast}\tilde{u}_{2}+\tilde{h}\qquad\text{in
}B_{1/2}^{4}$

for some harmonic 1-form $\tilde{h}$, where $\tilde{u}_{1}$ is a function and
$\tilde{u}_{2}$ is a $2$-form. Then $\tilde{u}_{1},\tilde{u}_{2}$ satisfy
equations

$\Delta^{2}\tilde{u}_{1}={\rm div}\tilde{K}+Af\qquad\text{in }B_{1/2}^{4},$
$\Delta\tilde{u}_{2}=dA\wedge du\qquad\text{in }B_{1/2}^{4}.$

To proceed, let us extend $V,w,\omega,F,A,B$ from $B_{1/2}^{4}$ to the whole
space $\mathbb{R}^{4}$ such that for some $C_{m}>0$ independent of $\epsilon$,
there hold

$\|V\|_{W^{1,2}(\mathbb{R}^{4})}+\|w\|_{L^{2}(\mathbb{R}^{4})}+\|\omega\|_{L^{2}(\mathbb{R}^{4})}+\|F\|_{L^{4/3,1}(\mathbb{R}^{4})}<C_{m}\epsilon,$
(2.1)
$\|A\|_{W^{2,2}(\mathbb{R}^{4})}+\|B\|_{W^{1,4/3}(\mathbb{R}^{4})}<C_{m}\epsilon$

(still denoted by the same notations). As to $f$, we simply let $f\equiv 0$ in
$\mathbb{R}^{4}\backslash B_{1/2}^{4}$ so that

$\|f\|_{M^{1,\lambda}(\mathbb{R}^{4})}\leq
C(\lambda)\|f\|_{M^{1,\lambda}(B_{1/2}^{4})}$

for some constant $C(\lambda)$ depending only on $\lambda$.

Put

$u_{11}=c\log\ast\left({\rm div}\tilde{K}\right),\qquad
u_{12}=c\log\ast(Af),\quad\quad u_{2}=I_{2}(dA\wedge du),$

where $c\log$ and $I_{2}=c|x|^{-2}$ denote the fundamental solutions of
$\Delta^{2}$ and $\Delta$ in $\mathbb{R}^{4}$ respectively, so that

$\Delta^{2}u_{11}={\rm
div}\tilde{K},\quad\quad\Delta^{2}u_{12}=Af,\quad\quad\Delta u_{2}=dA\wedge
du$

in $\mathbb{R}^{4}$. Then we find that $\Delta^{2}(u_{11}+u_{12}-u_{1})=0$ and
$\Delta(u_{2}-\tilde{u}_{2})=0$ in $B_{1/2}^{4}$ respectively. Thus, we obtain
the decomposition

(2.2) $Adu=du_{11}+du_{12}+d^{\ast}u_{2}+h$

for some biharmonic function $h$ in $B_{1/2}^{4}$.

Now we can use Riesz potential theory to estimate each term in equation (2.2).

For the first term $u_{11}$, note that

$|\nabla u_{11}|\lesssim I_{2}(|\tilde{K}|).$

$\tilde{K}$ contains 7 terms. For the first term $\nabla A\Delta u$, we have

$\left\|I_{2}(\nabla A\Delta
u)\right\|_{L^{4}(\mathbb{R}^{4})}\lesssim\|\nabla
A\|_{L^{4}(\mathbb{R}^{4})}\|\nabla^{2}u\|_{L^{2}(\mathbb{R}^{4})}\lesssim\epsilon\|\nabla^{2}u\|_{L^{2}(\mathbb{R}^{4})},$

where we used the boundedness of the operator
$I_{2}:L^{4/3}(\mathbb{R}^{4})\to L^{4}(\mathbb{R}^{4})$ and the estimate
(2.1). Similarly we can estimate the remaining six terms of $\nabla u_{11}$ to
find that

$\left\|\nabla
u_{11}\right\|_{L^{4}(\mathbb{R}^{4})}\lesssim\epsilon\left(\|\nabla^{2}u\|_{L^{2}(\mathbb{R}^{4})}+\|\nabla
u\|_{L^{4}(\mathbb{R}^{4})}\right).$

Since $u$ is extended to $\mathbb{R}^{4}$ in a bounded way from $B_{1/2}^{4}$,
we obtain

(2.3) $\left\|\nabla
u_{11}\right\|_{L^{4}(\mathbb{R}^{4})}\lesssim\epsilon\left(\|\nabla^{2}u\|_{L^{2}(B_{1/2}^{4})}+\|\nabla
u\|_{L^{4}(B_{1/2}^{4})}\right).$

To estimate $u_{12}$, we use the estimate

$|\nabla u_{12}|\lesssim I_{3}(|\nabla A||\nabla u|)$

and the boundedness of the operator (see Proposition A.1 in Appendix A)

$I_{3}:M^{1,\lambda}(\mathbb{R}^{4})\to
M_{\ast}^{\frac{4-\lambda}{1-\lambda},\lambda}(\mathbb{R}^{4})$

and find that $\nabla u_{12}\in
M_{\ast}^{\frac{4-\lambda}{1-\lambda},\lambda}(\mathbb{R}^{4})$ with estimate

$\left\|\nabla
u_{12}\right\|_{M_{\ast}^{\frac{4-\lambda}{1-\lambda},\lambda}(\mathbb{R}^{4})}\lesssim\|f\|_{M^{1,\lambda}(\mathbb{R}^{4})}\leq
C(\lambda)\|f\|_{M^{1,\lambda}(B_{1/2}^{4})}.$

In particular, for any $r>0$, we obtain by Hölder’s inequality that

$\|\nabla u_{12}\|_{L^{4}(B_{r}^{4})}\leq
C\left\|u_{12}\right\|_{L_{\ast}^{\frac{4-\lambda}{1-\lambda}}(B_{r}^{4})}r^{4\left(\frac{1}{4}-\frac{1-\lambda}{4-\lambda}\right)}\leq
C\|f\|_{M^{1,\lambda}(B_{1/2}^{4})}r^{\lambda}.$

Thus, $u_{12}\in M^{4,4\lambda}(\mathbb{R}^{4})$ with

(2.4) $\left\|\nabla u_{12}\right\|_{M^{4,4\lambda}(\mathbb{R}^{4})}\leq
C\|f\|_{M^{1,\lambda}(B_{1/2}^{4})}.$

The estimate of $u_{2}$ is totally the same as above, which gives

(2.5) $\left\|\nabla
u_{2}\right\|_{L^{4}(\mathbb{R}^{4})}\lesssim\epsilon\|\nabla
u\|_{L^{4}(\mathbb{R}^{4})}\lesssim\epsilon\left(\|\nabla^{2}u\|_{L^{2}(B_{1/2}^{4})}+\|\nabla
u\|_{L^{4}(B_{1/2}^{4})}\right).$

Since $h$ is biharmonic in $B_{1/2}^{4}$, for any $1/4>r>0$, we have

(2.6) $\|h\|_{L^{4}(B_{r}^{4})}\leq Cr\|h\|_{L^{4}(B_{1/2}^{4})}.$

Therefore, combining the above estimates (2.3), (2.4), (2.5) and (2.6) on
$u_{11}$, $u_{12}$, $u_{2}$ and $h$, we obtain, for any $1/4>r>0$,

$\displaystyle\|\nabla u\|_{L^{4}(B_{r}^{4})}$
$\displaystyle\leq\|h\|_{L^{4}(B_{r}^{4})}+\|\nabla
u_{11}\|_{L^{4}(B_{r}^{4})}+\|\nabla u_{12}\|_{L^{4}(B_{r}^{4})}+\|\nabla
u_{2}\|_{L^{4}(B_{r}^{4})}$ $\displaystyle\lesssim
r\|h\|_{L^{4}(B_{1/2}^{4})}+\epsilon\left(\|\nabla^{2}u\|_{L^{2}(B_{1/2}^{4})}+\|\nabla
u\|_{L^{4}(B_{1/2}^{4})}\right)+\|f\|_{M^{1,\lambda}(B_{1/2}^{4})}r^{\lambda}$
$\displaystyle\leq
C(r+\epsilon)\left(\|\nabla^{2}u\|_{L^{2}(B_{1/2}^{4})}+\|\nabla
u\|_{L^{4}(B_{1/2}^{4})}\right)+C\|f\|_{M^{1,\lambda}(B_{1/2}^{4})}r^{\lambda},$

where $C$ is a constant depending only on $m,\lambda$.

Similarly, repeating the above procedure, we can obtain the same estimate for
$\|\nabla^{2}u\|_{L^{2}(B_{r}^{4})}$. Hence, by writing

$U(r)\equiv\|\nabla
u\|_{L^{4}(B_{r}^{4})}+\|\nabla^{2}u\|_{L^{2}(B_{r}^{4})},$

we derive for any $0<r<1/4$

$U(r)\leq C(r+\epsilon)U(1/2)+C\|f\|_{M^{1,\lambda}(B_{1/2}^{4})}r^{\lambda}$

for some constant $C>0$ depending only on $m$ and $\lambda$. We first
$r_{0}\ll 1/4$ such that $2Cr_{0}\leq r_{0}^{(\lambda+1)/2}$, and then choose
$\epsilon\leq r_{0}$ so that

$U(r_{0})\leq
U(1/2)r_{0}^{(\lambda+1)/2}+C\|f\|_{M^{1,\lambda}(B_{1/2}^{4})}r_{0}^{\lambda}.$

Now using a standard scaling and iteration argument (see [10, Theorem 3.1] for
instance), we obtain, for any $k\geq 1$,

$U(r_{0}^{k})\leq
r_{0}^{\frac{\lambda+1}{2}}U(r_{0}^{k-1})+C\|f\|_{M^{1,\lambda}(B_{1/2}^{4})}r_{0}^{k\lambda},$

which implies that

$U(r)\leq C\left(U(1)+\|f\|_{M^{1,\lambda}(B_{1}^{4})}\right)r^{\lambda}$

for all $0<r<1/4$. This proves (1.9).

The Hölder continuity estimate (1.10) follows from a standard Morrey type
Dirichlet growth theorem, see Giaquinta [6] for example. This completes the
proof of Theorem 1.2.

##  3\. Proof of theorem 1.3

In this section we prove Theorem 1.3.

Suppose now $u\in W^{2,2}(B^{4})$ is a weak solution to equation (1.1) with
$f\in M^{p,p\lambda}(B_{1}^{4})$ for some $1<p<4/3$ and $0<\lambda<(4-3p)/p$.
Since the result is local and the problem is invariant under suitable scaling
(for details, see [10]), we can assume that

$\|V\|_{W^{1,2}(B_{1})}+\|w\|_{L^{2}(B_{1})}+\|\omega\|_{L^{2}(B_{1})}+\|F\|_{L^{4/3,1}(B_{1})}<\epsilon_{m}$

for some $\epsilon_{m}\ll 1$. Under this assumption, let us first improve the
result of Theorem 1.2 for later use.

1\. By Hölder’s inequality, we easily find that

(3.1) $M^{p,p\lambda}(B_{1}^{4})\subset M^{1,\lambda_{0}}(B_{1}^{4}),$

where

(3.2) $\lambda_{0}=\lambda+4(1-1/p).$

Thus by Theorem 1.2, we have

$\|\nabla
u\|_{M^{4,4\lambda_{0}}(B_{1/2}^{4})}+\|\nabla^{2}u\|_{M^{2,2\lambda_{0}}(B_{1/2}^{4})}\leq
C\left(\|u\|_{W^{2,2}(B^{4})}+\|f\|_{M^{1,\lambda_{0}}(B^{4})}\right)$

and $u\in C^{0,\lambda_{0}}(B_{1/2}^{4})$ with

(3.3) $\|u\|_{C^{0,\lambda_{0}}(B_{1/2}^{4})}\leq
C\left(\|u\|_{W^{2,2}(B^{4})}+\|f\|_{M^{1,\lambda_{0}}(B^{4})}\right).$

2\. Since $f\in L^{p}(B^{4})$ for $1<p<4/3$, by Guo-Xiang-Zheng [10], $u\in
W^{3,p_{1}}(B_{1/2}^{4})$ with

$\|u\|_{W^{3,p_{1}}(B_{1/2}^{4})}\leq
C\left(\|u\|_{W^{2,2}(B^{4})}+\|f\|_{L^{p}(B^{4})}\right),$

where

$p_{1}=4p/(4-p).$

3\. Due to the above $W^{3,p_{1}}$-regularity, we can then repeat the proof of
Theorem 1.2 to find that

$\|\nabla^{3}u\|_{M^{\frac{4}{3},\frac{4}{3}\lambda_{0}}(B_{1/2}^{4})}\leq
C\left(\|u\|_{W^{2,2}(B^{4})}+\|f\|_{M^{1,\lambda_{0}}(B^{4})}\right).$

Hence, we summarize the above result to conclude that

(3.4) $\|\nabla
u\|_{M^{4,4\lambda_{0}}(B_{1/2}^{4})}+\|\nabla^{2}u\|_{M^{2,2\lambda_{0}}(B_{1/2}^{4})}+\|\nabla^{3}u\|_{M^{\frac{4}{3},\frac{4}{3}\lambda_{0}}(B_{1/2}^{4})}\leq
CM,$

where

(3.5) $M=\|u\|_{W^{2,2}(B^{4})}+\|f\|_{M^{p,p\lambda}(B^{4})}.$

Now we can prove Theorem 1.3. Due to the second remark above, we know now that
$u\in W^{3,p_{1}}(B_{1/2}^{4})$. Our aim is to deduce its Morrey regularity.

Fix $x_{0}\in B_{1/4}^{4}$ and $0<r<1/4$. Split $A\Delta u=v+h$ in
$B_{r}(x_{0})$ with $v$ and $h$ satisfying

(3.6) $\begin{cases}\Delta v={\rm div}K+Af,&\text{in }B_{r}(x_{0}),\\\
v=0&\text{on }\partial B_{r}(x_{0}),\end{cases}$

and

(3.7) $\begin{cases}\Delta h=0&\text{in }B_{r}(x_{0}),\\\ h=A\Delta u&\text{on
}\partial B_{r}(x_{0}).\end{cases}$

1\. Estimate of $v$.

Write $p_{2}=4p/(4-2p)$ and $p_{3}=4p/(4-3p)$. By the duality of
$L^{p}$-space, we have

$\|v\|_{L^{p_{2}}(B_{r}(x_{0}))}=\sup_{\varphi\in{\mathcal{A}}}\int_{B_{r}(x_{0})}v\varphi,$

where

${\mathcal{A}}=\\{\varphi\in
L^{p_{2}^{\prime}}(B_{r}(x_{0});\mathbb{R}^{m}):\|\varphi\|_{p_{2}^{\prime}}\leq
1\\}$

and $p_{2}^{\prime}=p_{2}/(p_{2}-1)$ is the Hölder conjugate of $p_{2}$. For
any $\varphi\in{\mathcal{A}}$, the standard elliptic regularity theory implies
that there is a unique $\psi\in W^{2,p_{2}^{\prime}}\cap
W_{0}^{1,p_{2}^{\prime}}(B_{r}(x_{0}))$ satisfying $\Delta\psi=\varphi$ in
$B_{r}(x_{0})$, $\psi=0$ on $\partial B_{r}(x_{0})$, and moreover,

$\|\psi\|_{W^{2,p_{2}^{\prime}}(B_{r}(x_{0}))}\leq C(p).$

The Sobolev embedding theorem implies that

$W^{2,p_{2}^{\prime}}(B_{r}(x_{0}))\subset
W^{1,\frac{4p_{2}^{\prime}}{4-p_{2}^{\prime}}}(B_{r}(x_{0}))\subset
L^{p^{\prime}}(B_{r}(x_{0})).$

Hence

$\left\|\psi\right\|_{L^{p^{\prime}}(B_{r}(x_{0}))}+\left\|\nabla\psi\right\|_{L^{4p_{2}^{\prime}/(4-p_{2}^{\prime})}(B_{r}(x_{0}))}\leq
C(p).$

Now integrating by parts yields

$\int_{B_{r}(x_{0})}v\varphi=\int_{B_{r}(x_{0})}v\Delta\psi=-\int_{B_{r}(x_{0})}K\nabla\psi+\int_{B_{r}(x_{0})}Af\psi=:I_{1}+I_{2}.$

Using Hölder’s inequality we easily obtain

$|I_{2}|\leq
C\|f\|_{L^{p}(B_{r}(x_{0}))}\|\psi\|_{L^{p^{\prime}}(B_{r}(x_{0}))}\leq
C\|f\|_{M^{p,p\lambda}(B^{4})}r^{\lambda}.$

Here we used the assumption that $f\in M^{p,p\lambda}(B^{4})$ in the last
inequality. To estimate $I_{1}$, we use Hölder’s inequality again to infer
that

$\int_{B_{r}(x_{0})}|\nabla A\Delta u\nabla\psi|\leq\|\nabla A\|_{4}\|\|\Delta
u\|_{L^{p_{2}}(B_{r}(x_{0}))}\left\|\nabla\psi\right\|_{L^{4p_{2}^{\prime}/(4-p_{2}^{\prime})}(B_{r}(x_{0}))}\leq
C\epsilon\|\Delta u\|_{L^{p_{2}}(B_{r}(x_{0}))}.$

Other terms of $K$ can be estimated similarly. This gives

$|I_{1}|\leq C\epsilon\left(\|\Delta u\|_{L^{p_{2}}(B_{r}(x_{0}))}+\|\nabla
u\|_{L^{p_{3}}(B_{r}(x_{0}))}\right)$

Hence combining the above two estimates on $I_{1}$ and $I_{2}$ gives

$\|v\|_{L^{p_{2}}(B_{r}(x_{0}))}\leq C\epsilon\left(\|\Delta
u\|_{L^{p_{2}}(B_{r}(x_{0}))}+\|\nabla
u\|_{L^{p_{3}}(B_{r}(x_{0}))}\right)+C\|f\|_{M^{p,p\lambda}(B^{4})}r^{\lambda}.$

We can estimate the norm of $\nabla v$ similarly. Note that

(3.8) $\|\nabla
v\|_{L^{p_{1}}(B_{r}(x_{0}))}=\sup_{\varphi\in{\mathcal{B}}}\int_{B_{r}(x_{0})}\nabla
v\cdot\varphi,$

where

${\mathcal{B}}=\left\\{\varphi\in
L^{p_{1}^{\prime}}(B_{r}(x_{0}),\mathbb{R}^{m}\otimes\mathbb{R}^{4}):\|\varphi\|_{L^{p_{1}^{\prime}}(B_{r}(x_{0}))}\leq
1\right\\}.$

Using Helmholtz decomposition, for each $\varphi\in{\mathcal{B}}$, there exist
$\psi\in W_{0}^{1,p_{1}^{\prime}}(B_{r}(x_{0}),\mathbb{R}^{m})$ and $\xi\in
L^{p_{1}^{\prime}}(B_{r}(x_{0}),\mathbb{R}^{m}\otimes\mathbb{R}^{4})$ such
that $\varphi=\nabla\psi+\xi$ and ${\rm div}\xi=0$ in $B_{r}(x_{0})$.
Substituting this decomposition into (3.8) gives

$\int_{B_{r}(x_{0})}\nabla v\cdot\varphi=\int_{B_{r}(x_{0})}\nabla
v\cdot\nabla\psi$

in view of the boundary condition $v=0$ on $\partial B_{r}(x_{0})$. Working as
in the above for $\|v\|_{L^{p_{2}}(B_{r}(x_{0}))}$, we deduce

$\|\nabla v\|_{L^{p_{1}}(B_{r}(x_{0}))}\leq C\epsilon\left(\|\Delta
u\|_{L^{p_{2}}(B_{r}(x_{0}))}+\|\nabla
u\|_{L^{p_{3}}(B_{r}(x_{0}))}\right)+C\|f\|_{M^{p,p\lambda}(B^{4})}r^{\lambda}.$

Therefore, we conclude that

(3.9) $\displaystyle\qquad\|\nabla
v\|_{L^{p_{1}}(B_{r}(x_{0}))}+\|v\|_{L^{p_{2}}(B_{r}(x_{0}))}$
$\displaystyle\leq C\epsilon\left(\|\Delta
u\|_{L^{p_{2}}(B_{r}(x_{0}))}+\|\nabla
u\|_{L^{p_{3}}(B_{r}(x_{0}))}\right)+C\|f\|_{M^{p,p\lambda}(B^{4})}r^{\lambda}.$

2\. Estimate of $h$.

Since $A\Delta u\in W^{1,4/3}(B_{r}(x_{0}))$, using a standard result (see
e.g. Lemma C.1 of [10]) we obtain that $h\in W^{1,4/3}(B_{r}(x_{0}))$ with

$\|h\|_{W^{1,4/3}(B_{r}(x_{0}))}\leq C\|A\Delta
u\|_{W^{1,4/3}(B_{r}(x_{0}))}\leq C\Phi(r)r^{\lambda_{0}},$

where $\lambda_{0}$ are defined as in (3.2) and

$\Phi(r)=\|\nabla\Delta u\|_{L^{p_{1}}(B_{r}(x_{0}))}+\|\Delta
u\|_{L^{p_{2}}(B_{r}(x_{0}))}+\|\nabla u\|_{L^{p_{3}}(B_{r}(x_{0}))}.$

In particular, this implies that

$\|h\|_{W^{1,4/3}(B_{r}(x_{0}))}\leq CMr^{\lambda_{0}},$

where $M$ is defined as in (3.5). Since $h$ is harmonic in $B_{r}(x_{0})$,
this implies the interior growth estimate

$\|h\|_{L^{p_{2}}(B_{r/2}(x_{0}))}\leq
Cr^{\frac{4}{p_{2}}}\|h\|_{L^{\infty}(B_{r/2}(x_{0}))}\leq
Cr^{\frac{4}{p_{2}}}\left(\fint_{B_{r}(x_{0})}|h|^{2}\right)^{\frac{1}{2}}\leq
CMr^{\lambda}$

and similarly

$\|\nabla h\|_{L^{p_{1}}(B_{r/2}(x_{0}))}\leq CMr^{\lambda}.$

3\. before conclusion, we need to estimate $\nabla u$. Applying the Sobolev
embedding $W^{2,p_{2}}(B_{r/2}^{4})\subset W^{1,p_{3}}(B_{r/2}^{4})$ we find
that

$\displaystyle\|\nabla u\|_{L^{p_{3}}(B_{r/4}(x_{0}))}$ $\displaystyle\leq
C\|\nabla^{2}u\|_{L^{p_{2}}(B_{r/4}(x_{0}))}+Cr^{-2}\|u-u_{B_{r/4}(x_{0})}\|_{L^{p_{2}}(B_{r/4}(x_{0}))}$
$\displaystyle\leq C\|\Delta
u\|_{L^{p_{2}}(B_{r/2}(x_{0}))}+Cr^{-4(1-1/p)}\|u-u_{B_{r/2}(x_{0})}\|_{L^{\infty}(B_{r/2}(x_{0}))}$

for some constant $C>0$ independent of $r$. Then use the Hölder estimate (3.3)
we get

(3.10) $\|\nabla u\|_{L^{p_{3}}(B_{r/4}(x_{0}))}\leq
C\left(\|v\|_{L^{p_{2}}(B_{r/2}(x_{0}))}+\|h\|_{L^{p_{2}}(B_{r/2}(x_{0}))}\right)+CMr^{\lambda}.$

Finally, by a slight refinement of the above argument, combining the estimates
(3.9) of $v$ and $h$ and the estimate (3.10) of $\nabla u$, we conclude that

$\Phi(r/2)\leq C\epsilon\Phi(r)+CMr^{\lambda}.$

At this moment we use Simon’s iteration lemma (see e.g. Lemma A.7 of [17]) to
infer that there are $\epsilon_{0}$ and $r_{0}(\lambda,p)>0$ sufficiently
small such that for all $r\leq r_{0}$, we have

$\Phi(r)\leq CMr^{\lambda}.$

This completes the proof.

##  Appendix A Morrey spaces and Riesz potential

Let $\Omega\subset\mathbb{R}^{n}$ be a domain and let $1\leq p<\infty$ and
$0\leq s\leq n$. Denote by $B_{r}(x)$ the ball in $\mathbb{R}^{n}$ centered at
$x\in\mathbb{R}^{n}$ with radius $r$. The Morrey space $M^{p,s}(\Omega)$
consists of functions $f\in L^{p}(\Omega)$ such that

$\|f\|_{M^{p,s}(\Omega)}\equiv\sup_{x\in\Omega,0<r<\text{\rm
diam}\Omega}r^{-s/p}\|f\|_{L^{p}(B_{r}(x)\cap\Omega)}<\infty.$

The weak Morrey space $M_{\ast}^{p,s}(\Omega)$ consists of weakly
$p$-integrable functions $f\in L_{\ast}^{p}(\Omega)$ such that

$\|f\|_{M_{\ast}^{p,s}(\Omega)}\equiv\sup_{x\in D,0<r<\text{\rm
diam}\Omega}r^{-s/p}\|f\|_{L_{\ast}^{p}(B_{r}(x)\cap\Omega)}<\infty.$

Up to a constant multiplier, the upper bound of $r$ in the supremum above can
be also replaced by $(\text{\rm diam}\Omega)/2$ for instance. For more
details, see e.g. [2]. The following potential theory is due to Adams [1].

######  Proposition A.1.

Let $0<\alpha<n$ and $I_{\alpha}=|x|^{\alpha-n}$. Let $1\leq p<n/\alpha$ and
$0\leq\lambda<n$. Then the following linear operators are bounded:

(1) for $1<p<(n-\lambda)/\alpha$,

$I_{\alpha}:M^{p,\lambda}(\mathbb{R}^{n})\to
M^{\frac{(n-\lambda)p}{n-\lambda-\alpha p},\lambda}(\mathbb{R}^{n}),$

(2) for $0<\alpha<n-\lambda$,

$I_{\alpha}:M^{1,\lambda}(\mathbb{R}^{n})\to
M_{\ast}^{\frac{n-\lambda}{n-\lambda-\alpha},\lambda}(\mathbb{R}^{n}).$

## References

  * [1] D.R. Adams, _A note on Riesz potentials_. Duke Math. J. 42 (1975), no. 4, 765-778.
  * [2] R.C. Adams and J.F. Fournier, _Sobolev spaces._ Second edition. Pure and Applied Mathematics (Amsterdam), 140. Elsevier/Academic Press, Amsterdam, 2003.
  * [3] S.-Y. A. Chang, L. Wang and P.C. Yang, _A regularity theory of biharmonic maps._ Commun. Pure Appl. Math. 52(9) (1999), 1113-1137.
  * [4] H. Du, Y. Kang and J. Wang, _Morrey regularity theory of Rivière’s equation._ To appear in Proc. Amer. Math. Soc. 2022.
  * [5] A. Gastel, _The extrinsic polyharmonic map heat flow in the critical dimension._ Adv. Geom. 6 (2006), no. 4, 501-521.
  * [6] M. Giaquinta, _Multiple integrals in the calculus of variations and nonlinear elliptic systems._ Annals of Mathematics Studies, 105. Princeton University Press, Princeton, NJ, 1983.
  * [7] C.-Y. Guo, C. Wang and C.-L. Xiang, _$L^{p}$ -regularity for fourth order elliptic systems with antisymmetric potentials in higher dimensions._ Calc. Var. Partial Differential Equations 62 (2023), no. 1, Paper No. 31.
  * [8] C.-Y. Guo and C.-L. Xiang, _Regularity of solutions for a fourth order linear system via conservation law_. J. Lond. Math. Soc. (2) 101 (2020), no. 3, 907-922.
  * [9] C.-Y. Guo and C.-L. Xiang, _Regularity of weak solutions to higher order elliptic systems in critical dimensions_. Tran. Amer. Math. Soc., 374 (2021), no. 5, 3579-3602.
  * [10] C.-Y. Guo, C.-L. Xiang and G.-F. Zheng, _The Lamm-Riviere system I: $L^{p}$ regularity theory_. Calc. Var. Partial Differential Equations 60 (2021), no. 4, Paper No. 125, 23 pp.
  * [11] C.-Y. Guo, C.-L. Xiang and G.-F. Zheng, _$L^{p}$ regularity theory for even order elliptic systems with antisymmetric first order potentials._ J. Math. Pures Appl. (9) 165 (2022), 286-324.
  * [12] T. Lamm, _Heat flow for extrinsic biharmonic maps with small initial energy._ Ann. Global Anal. Geom. 26 (2004), no. 4, 369-384.
  * [13] T. Lamm and T. Rivière, _Conservation laws for fourth order systems in four dimensions._ Comm. Partial Differential Equations 33 (2008), 245-262.
  * [14] P. Laurain and T. Rivière, _Angular energy quantization for linear elliptic systems with antisymmetric potentials and applications._ Anal. PDE 7 (2014), no. 1, 1-41.
  * [15] T. Rivière, _Conservation laws for conformally invariant variational problems._ Invent. Math. 168 (2007), 1-22.
  * [16] T. Rivière and M. Struwe, _Partial regularity for harmonic maps and related problems._ Comm. Pure Appl. Math. 61 (2008), 451-463.
  * [17] B. Sharp and P. Topping, _Decay estimates for Rivière’s equation, with applications to regularity and compactness_. Trans. Amer. Math. Soc. 365 (2013), no. 5, 2317-2339.
  * [18] M. Struwe, _Partial regularity for biharmonic maps, revisited._ Calc. Var. Partial Differential Equations 33 (2008), 249-262.
  * [19] C. Wang, _Heat flow of biharmonic maps in dimensions four and its application._ Pure Appl. Math. Q. 3 (2007), no. 2, Special Issue: In honor of Leon Simon. Part 1, 595-613.
  * [20] C. Wang, _Remarks on approximate harmonic maps in dimension two._ Calc. Var. Partial Differential Equations 56 (2017), no. 2, Paper No. 23, 24 pp.
  * [21] C.Y. Wang, _Biharmonic maps from $R^{4}$ into a Riemannian manifold._ Math. Z. 247 (2004), 65-87.
  * [22] C.Y. Wang, _Stationary biharmonic maps from Rm into a Riemannian manifold._ Comm. Pure Appl. Math. 57 (2004), 419-444.