# 1+2 dimensional radially symmetric wave maps revisit.

Yi Zhou School of Mathematics Science, Fudan University, Shanghai, P. R. China
(yizhou@fudan.edu.cn).

###### Abstract

The author gives an alternative and simple proof of the global existence of
smooth solutions to the Cauchy problem for wave maps from the 1+2-dimensional
Minkowski space to an arbitrary compact smooth Riemannian manifold without
boundary, for arbitrary smooth, radially symmetric data. the author can also
treat non-compact manifold under some additional assumptions which generalize
the existing ones.

Keywords: Cauchy problem, wave maps, global smooth solution

##  1 Main result

Let $N$ be a smooth Riemannian k-manifold without boundary. With no loss of
generality, we assume that $N\subset\mathbb{R}^{n}$, isometrically. We
consider wave maps
$\Phi=(\Phi^{1},\cdots,\Phi^{n})=\Phi(t,x):\mathbb{R}\times\mathbb{R}^{2}\rightarrow
N\in\mathbb{R}^{n}$ , satisfying the equation

$\Box\Phi=\Phi_{tt}-\Delta\Phi=B(\Phi)(\partial_{\alpha}\Phi,\partial^{\alpha}\Phi)\perp
T_{\Phi}N,$ (1.1)

where $B$ denotes the second fundamental form of $N$. Writing
$z=(t,x)=(x^{\alpha})_{0\leq\alpha\leq 2}$, we also let
$\partial_{\alpha}=\frac{\partial}{\partial x^{\alpha}},\ \alpha=0,1,2.$ We
raise and lower indices with the Minkowski metric
$\eta=(\eta_{\alpha\beta})=(\eta^{\alpha\beta})=diag(-1,1,1)$ and tacitly sum
over repeated indices.

Due to its mathematical difficulty and important physical background, the
topic of wave maps has experienced an incredible advancement in the past
several decades. It has at least two physical motivations to study wave maps.
One is the nonlinear $\sigma$-model which deals with the case $N$ is a sphere
and on the other hand, vacuum Einstein equations with $U(1)\times\mathbb{R}$
symmetries reduce to a radially symmetric wave maps from 1+2 dimensional
Minkowski space to the hyperbolic plane. It is Prof. Gu [3] who first gave the
regularity result in $1+1$ dimensional case. For an up to date account of the
full developments, we refer to the monograph of Geba, D.A. and Grillakis M.G.
[2], as well as [4], [5], and [8].

The main purpose of this paper is to give an alternative proof of the global
existence for 1+2 dimensional wave maps with radial symmetry. We may assume
that TN is parallelizable. Let $\bar{e}_{1},\cdots,\bar{e}_{k}$ be a smooth
orthonormal frame field such that at any point $p\in N$ the vectors
$\bar{e}_{1}(p),\cdots,\bar{e}_{k}(p)$ form an orthonormal basis for
$T_{p}N$.We make the following assumption on the target manifold.

(H1): $\nabla_{p}\bar{e}_{i}(p)$ is uniformly bounded,

(H2): $\nabla_{p}B(p)$ is uniformly bounded.

We note that those assumptions are in particular satisfied for compact
manifold.

The main result of this paper is the following.

######  Theorem 1.1.

Let $N\subset\mathbb{R}^{n}$ be a smooth Riemannian manifold without boundary.
Assume (H1) (H2) satisfied for $N$. Then for any radially symmetric data

$(\Phi_{0},\Phi_{1})=(\Phi_{0}(r),\Phi_{1}(r))\in
C^{\infty}(\mathbb{R}^{2},TN),\ r=|x|,$ (1.2)

there exists a unique, smooth solution $\Phi=\Phi(t,r)$ to the Cauchy problem
(1.1)(1.2), defined for all time.

Our result slightly generalizes the work of Christodoulou-Tahvildar-Zadeh [1]
and Struwe [6],[7].

As usual, the proof of Theorem1.1 is divided into two main steps. The first
step is to show that small energy implies regularity and the second step is to
show that energy can not concentrate. In this paper, we give an alternative
proof of the first step which is totally different from the work of
Christodoulou-Tahvildar-Zadeh [1] and the proof of the second step can be
found in Struwe [6],[7]. In the work of Christodoulou-Tahvildar-Zadeh [1], the
first step is achieved by a H$\ddot{o}$lder estimate using the fundamental
solution of the 2-dimensional radially symmetric wave operator which is quite
complicated, While we rely on an energetic argument which we call a new div
curl lemma and has potential to work for quasilinear problems. This will be
pursued in our future work.

##  2 Intrinsic setting

Let $\Phi=\Phi(t,r)$ be a smooth radially symmetric wave maps, then (1.1)(1.2)
can be rewriting as

$\displaystyle\Phi_{tt}-\Phi_{rr}-\frac{\Phi_{r}}{r}=4B(\Phi)(\Phi_{u},\Phi_{v}),\
(t>0)$ (2.1)

$\displaystyle\Phi=\Phi_{0}(r),\ \Phi_{t}=\Phi_{1}(r),\ (t=0)$ (2.2)

we denote $u=t-r,\ v=t+r,$ then

$\displaystyle
t=\frac{1}{2}(u+v),\partial_{t}=\partial_{u}+\partial_{v},\partial_{u}=\frac{1}{2}(\partial_{t}-\partial_{r}),$
$\displaystyle
r=\frac{1}{2}(v-u),\partial_{r}=\partial_{v}-\partial_{u},\partial_{v}=\frac{1}{2}(\partial_{t}+\partial_{r}),$
$\displaystyle\Phi_{u}=\partial_{u}\Phi,\Phi_{v}=\partial_{v}\Phi,\Phi_{tt}=\partial_{t}^{2}\Phi,\
etc.$

and we have the energy conservation law

$\displaystyle
E(t)=\frac{1}{2}\int_{0}^{\infty}r(|\Phi_{t}|^{2}+|\Phi_{r}|^{2})dr=E_{0}=\frac{1}{2}\int_{0}^{\infty}r(|\Phi_{0}^{\prime}|^{2}+|\Phi_{1}|^{2})dr.$
(2.3)

Let $D$ be the pull-back covariant derivative in $u^{*}TN$, we may write
equation(2.1) as

$\displaystyle D_{t}\Phi_{t}-\frac{1}{r}D_{r}(r\Phi_{r})=0.$ (2.4)

From $(\bar{e}_{i})_{1\leq i\leq k}$, we obtain a frame
$e_{i}=R_{ij}(\bar{e}_{j}\circ\Phi),1\leq i\leq k$ for the pull-back bundle,
where $R=R(z)=(R_{ij})$ may be any smooth map. Denoting

$\displaystyle D_{\alpha}e_{i}=A_{i\alpha}^{j}e_{j},\ 0\leq\alpha\leq 2$ (2.5)

with a matrix-valued connection 1-form $A=A_{\alpha}dx^{\alpha}$. We compute
the curvature $F$ of $D$ via the commutation relation

$\displaystyle
D_{\alpha}D_{\beta}e_{i}-D_{\beta}D_{\alpha}e_{i}=D_{\alpha}(A_{i\beta}^{j}e_{j})-D_{\beta}(A_{i\alpha}^{j}e_{i})$
$\displaystyle=(\partial_{\alpha}A_{i\beta}^{k}-\partial_{\beta}A_{i\alpha}^{k}-A_{j\alpha}^{k}A_{i\beta}^{j}-A_{j\beta}^{k}A_{i\alpha}^{j})e_{k}$
$\displaystyle=F_{i\alpha\beta}^{k}e_{k}$

or, more precisely

$\displaystyle dA+\frac{1}{2}[A,A]=F.$ (2.6)

In the radially symmetric case, we may choose $R=R(t,r),A=A_{0}dt+A_{1}dr$.
Following Struwe [7], we may impose ”exponential gauge” condition $A_{1}=0$,
which yields the relation

$\displaystyle*dA=-\partial_{r}A_{0}=F_{01}.$ (2.7)

If we normalize $A_{0}(t,\infty)=0,\forall t$, from this relation we obtain

$\displaystyle A_{0}=\int_{r}^{+\infty}F_{01}ds.$ (2.8)

By (H1), we get

$\displaystyle|F_{01}|\leq C|d\Phi|^{2}.$ (2.9)

Thus we deduce the estimate

$\displaystyle|A_{0}|\lesssim\int_{r}^{+\infty}|d\Phi|^{2}ds\lesssim
E_{0}r^{-1}.$ (2.10)

Let

$\displaystyle\Phi_{t}=q_{0}^{i}e_{i},\ \Phi_{r}=q_{1}^{i}e_{i}.$ (2.11)

Using the notation

$\displaystyle D_{\alpha}\partial_{\beta}\Phi$
$\displaystyle=D_{\alpha}(q_{\beta}^{i}e_{i})=(\partial_{\alpha}q_{\beta}^{j}+A_{i\beta}^{j}q_{\beta}^{i})e_{j}$
$\displaystyle=(D_{\alpha}q_{\beta}^{i})^{j}e_{j},$ (2.12)

we then may write equation (2.4) in the form

$\displaystyle
D_{t}q_{0}-\frac{1}{r}D_{r}(rq_{1})=\partial_{t}q_{0}+A_{0}q_{0}-\frac{1}{r}\partial_{r}(rq_{1})=0.$
(2.13)

Moreover, we have the commutation relation

$\displaystyle D_{r}q_{0}=D_{t}q_{1}=\partial_{t}q_{1}+A_{0}q_{1}.$ (2.14)

Taking inner product in $\mathbb{R}^{n}$ of (2.13) with $-r^{1-\alpha}q_{1}$,
we get

$\displaystyle 0$
$\displaystyle=-r^{1-\alpha}q_{1}\partial_{t}q_{0}-r^{1-\alpha}q_{1}(A_{0}q_{0})+r^{1-\alpha}q_{1}\partial_{r}q_{1}+r^{-\alpha}q_{1}^{2}$
$\displaystyle=-r^{1-\alpha}q_{1}\partial_{t}q_{0}-r^{1-\alpha}q_{1}(A_{0}q_{0})+\left(r^{1-\alpha}q_{1}\partial_{r}q_{1}+\frac{1}{2}(1-\alpha)r^{-\alpha}q_{1}^{2}\right)+\frac{1}{2}(1+\alpha)r^{-\alpha}q_{1}^{2}$
$\displaystyle=-r^{1-\alpha}q_{1}\partial_{t}q_{0}-r^{1-\alpha}q_{1}(A_{0}q_{0})+\partial_{r}\left(\frac{1}{2}r^{1-\alpha}q_{1}^{2}\right)+\frac{1}{2}(1+\alpha)r^{-\alpha}q_{1}^{2}.$
(2.15)

Taking inner product in $\mathbb{R}^{n}$ of (2.14) with $-r^{1-\alpha}q_{0}$,
we get

$\displaystyle 0$
$\displaystyle=-r^{1-\alpha}q_{0}\partial_{t}q_{1}-r^{1-\alpha}q_{0}(A_{0}q_{1})+r^{1-\alpha}q_{0}\partial_{r}q_{0}$
$\displaystyle=-r^{1-\alpha}q_{0}\partial_{t}q_{1}-r^{1-\alpha}q_{0}(A_{0}q_{1})+\left(r^{1-\alpha}q_{0}\partial_{r}q_{0}+\frac{1}{2}(1-\alpha)r^{-\alpha}q_{0}^{2}\right)-\frac{1}{2}(1-\alpha)r^{-\alpha}q_{0}^{2}$
$\displaystyle=-r^{1-\alpha}q_{0}\partial_{t}q_{1}-r^{1-\alpha}q_{0}(A_{0}q_{1})+\partial_{r}\left(\frac{1}{2}r^{1-\alpha}q_{0}^{2}\right)-\frac{1}{2}(1-\alpha)r^{-\alpha}q_{0}^{2}.$
(2.16)

Combining (2) (2) and the antisymmetry of $A_{0}$, we know

$\displaystyle-\partial_{t}\left(r^{1-\alpha}q_{0}q_{1}\right)+\partial_{r}\left(\frac{1}{2}r^{1-\alpha}(q_{0}^{2}+q_{1}^{2})\right)-\frac{1}{2}(1-\alpha)r^{-\alpha}q_{0}^{2}+\frac{1}{2}(1+\alpha)r^{-\alpha}q_{1}^{2}$
$\displaystyle=r^{1-\alpha}q_{1}(A_{0}q_{0})+r^{1-\alpha}q_{0}(A_{0}q_{1})=0.$
(2.17)

Define $Q_{0}(r)=\int_{0}^{r}\xi^{-\sigma}q_{1}(\xi)d\xi$,
$Q(r)=r^{-\alpha+\sigma}Q_{0}$, $\sigma<\alpha$. And by (2.14), we have

$\displaystyle Q\partial_{t}q_{0}$
$\displaystyle=\partial_{t}(Qq_{0})-\partial_{t}Qq_{0}=\partial_{t}(r^{-\alpha+\sigma}Q_{0}q_{0})-r^{-\alpha+\sigma}\partial_{t}Q_{0}q_{0}$
$\displaystyle=\partial_{t}(r^{-\alpha+\sigma}Q_{0}q_{0})-\frac{q_{0}}{r^{\alpha-\sigma}}\int_{0}^{r}\xi^{-\sigma}\partial_{t}q_{1}(\xi)d\xi$
$\displaystyle=\partial_{t}(r^{-\alpha+\sigma}Q_{0}q_{0})-\frac{q_{0}}{r^{\alpha-\sigma}}\int_{0}^{r}\xi^{-\sigma}[\partial_{\xi}q_{0}(\xi)-A_{0}q_{1}]d\xi$
$\displaystyle=\partial_{t}(r^{-\alpha+\sigma}Q_{0}q_{0})+\frac{q_{0}}{r^{\alpha-\sigma}}\int_{0}^{r}\xi^{-\sigma}A_{0}q_{1}d\xi-\frac{q_{0}(r)}{r^{\alpha-\sigma}}\left[\left.\xi^{-\sigma}q_{0}(\xi)\right|_{0}^{r}+\sigma\int_{0}^{r}\xi^{-\sigma-1}q_{0}(\xi)d\xi\right]$
$\displaystyle=\partial_{t}(r^{-\alpha+\sigma}Q_{0}q_{0})+\frac{q_{0}}{r^{\alpha-\sigma}}\int_{0}^{r}\xi^{-\sigma}A_{0}q_{1}d\xi-\frac{q_{0}^{2}}{r^{\alpha}}-\frac{\sigma
q_{0}}{r^{\alpha-\sigma}}\int_{0}^{r}\xi^{-\sigma-1}q_{0}(\xi)d\xi$ (2.18)

$\displaystyle-\partial_{r}q_{1}Q$
$\displaystyle=-\partial_{r}(q_{1}Q)+q_{1}\partial_{r}Q=-\partial_{r}(q_{1}Q)+q_{1}\partial_{r}(r^{-\alpha+\sigma}Q_{0})$
$\displaystyle=-\partial_{r}(q_{1}Q)+(-\alpha+\sigma)r^{-\alpha+\sigma-1}q_{1}Q_{0}+r^{-\alpha+\sigma}q_{1}\partial_{r}Q_{0}$
$\displaystyle=-\partial_{r}(q_{1}r^{-\alpha+\sigma}Q_{0})+(-\alpha+\sigma)\frac{q_{1}Q}{r}+\dfrac{q_{1}^{2}}{r^{\alpha}}$
(2.19) $\displaystyle\frac{q_{1}}{r}Q$
$\displaystyle=\frac{q_{1}}{r^{\alpha+1-\sigma}}Q_{0}=\frac{1}{r^{\alpha+1-2\sigma}}\partial_{r}Q_{0}Q_{0}$
$\displaystyle=\frac{1}{2}\partial_{r}(r^{2\sigma-1-\alpha}Q_{0}^{2})-\frac{2\sigma-1-\alpha}{2}r^{2\sigma-2-\alpha}Q_{0}^{2}$
(2.20)

Combining (2) (2), we can get

$\displaystyle-\partial_{r}q_{1}Q-\frac{q_{1}}{r}Q$
$\displaystyle=-\partial_{r}(q_{1}r^{-\alpha+\sigma}Q_{0})+(-\alpha+\sigma-1)\frac{q_{1}Q}{r}+\dfrac{q_{1}^{2}}{r^{\alpha}}$
$\displaystyle=-\partial_{r}(q_{1}r^{-\alpha+\sigma}Q_{0})+\dfrac{q_{1}^{2}}{r^{\alpha}}$
$\displaystyle+\frac{-\alpha+\sigma-1}{2}\partial_{r}(r^{2\sigma-1-\alpha}Q_{0}^{2})-\frac{(2\sigma-1-\alpha)(-\alpha+\sigma-1)}{2}r^{2\sigma-2-\alpha}Q_{0}^{2}$
(2.21)

From (2) and (2), taking inner product in $\mathbb{R}^{n}$ of (2.13) with
$-\dfrac{Q}{2}$, we know

$\displaystyle-\frac{1}{2}$
$\displaystyle\partial_{t}(r^{-\alpha+\sigma}Q_{0}q_{0})+\frac{1}{2}\frac{q_{0}^{2}}{r^{\alpha}}-\frac{1}{2}\dfrac{q_{1}^{2}}{r^{\alpha}}+\frac{1}{2}\partial_{r}(q_{1}r^{-\alpha+\sigma}Q_{0})$
$\displaystyle+\frac{\alpha-\sigma+1}{4}\partial_{r}(r^{2\sigma-1-\alpha}Q_{0}^{2})+\frac{(1+\alpha-2\sigma)(\alpha-\sigma+1)}{4}r^{2\sigma-2-\alpha}Q_{0}^{2}$
$\displaystyle=\frac{q_{0}}{2r^{\alpha-\sigma}}\int_{0}^{r}\xi^{-\sigma}A_{0}q_{1}d\xi-\frac{\sigma
q_{0}}{2r^{\alpha-\sigma}}\int_{0}^{r}\xi^{-\sigma-1}q_{0}(\xi)d\xi+\frac{1}{2}r^{-\alpha+\sigma}A_{0}q_{0}Q_{0}$
(2.22)

By (2) and (2) , taking inner product in $\mathbb{R}^{n}$ of (2.13) with
$-(\dfrac{Q}{2}+r^{1-\alpha}q_{1})$, we obtain

$\displaystyle-\partial_{t}\left(r^{1-\alpha}q_{0}q_{1}\right)+\partial_{r}\left(\frac{1}{2}r^{1-\alpha}(q_{0}^{2}+q_{1}^{2})\right)+\frac{\alpha}{2}\frac{q_{0}^{2}+q_{1}^{2}}{r^{\alpha}}$
$\displaystyle-\frac{1}{2}\partial_{t}(r^{-\alpha+\sigma}Q_{0}q_{0})+\frac{1}{2}\partial_{r}(q_{1}r^{-\alpha+\sigma}Q_{0})$
$\displaystyle+\frac{\alpha-\sigma+1}{4}\partial_{r}(r^{2\sigma-1-\alpha}Q_{0}^{2})+\frac{(1+\alpha-2\sigma)(\alpha-\sigma+1)}{4}r^{2\sigma-2-\alpha}Q_{0}^{2}$
$\displaystyle=\frac{q_{0}}{2r^{\alpha-\sigma}}\int_{0}^{r}\xi^{-\sigma}A_{0}q_{1}d\xi-\frac{\sigma
q_{0}}{2r^{\alpha-\sigma}}\int_{0}^{r}\xi^{-\sigma-1}q_{0}(\xi)d\xi+\frac{1}{2}r^{-\alpha+\sigma}A_{0}q_{0}Q_{0}.$
(2.23)

##  3 New div curl Lemma

The purpose of this section is to prove the following lemma3.1 and lemma3.2.
We call Lemma3.2 a new div curl Lemma.

######  Lemma 3.1.

Suppose that

$\begin{cases}\partial_{u}F^{11}+\partial_{v}F^{12}=G^{1}\ \ (r\geq 0,0\leq
t\leq T)\\\ \partial_{u}F^{21}-\partial_{v}F^{22}=G^{2}\ \ (r\geq 0,0\leq
t\leq T),\end{cases}$ (3.1)

where $F^{11},F^{12},F^{21},F^{22}$ are all nonnegative. Moreover,

$\displaystyle r=0:\ F^{11}-F^{12}=F^{21}+F^{22}=0.$ (3.2)

Then there hold

$\displaystyle\sup_{u}\int_{|u|}^{2T-u}F^{11}(u,v)du+\sup_{v}\int_{-v}^{min\\{v,2T-v\\}}F^{12}(u,v)du$
$\displaystyle\lesssim\int_{0}^{\infty}(F^{11}+F^{12})(0,r)dr+\int_{0}^{T}\int_{0}^{\infty}|G_{1}|(t,r)drdt,$
(3.3)
$\displaystyle\sup_{u}\int_{|u|}^{2T-u}F^{21}(u,v)du+\sup_{v}\int_{-v}^{min\\{v,2T-v\\}}F^{22}(u,v)du$
$\displaystyle\lesssim\int_{0}^{\infty}(F^{21}+F^{22})(T,r)dr+\int_{0}^{\infty}(F^{21}+F^{22})(0,r)dr+\int_{0}^{T}\int_{0}^{\infty}|G_{2}|(t,r)drdt.$

Proof We only prove

$\displaystyle\sup_{u}\int_{|u|}^{2T-u}F^{11}(u,v)du$
$\displaystyle\lesssim\int_{0}^{\infty}(F^{11}+F^{12})(0,r)dr+\int_{0}^{T}\int_{0}^{\infty}|G_{1}|(t,r)drdt.$

The other estimates are similar. Noting(3.2), we have

$\displaystyle\partial_{u}\int_{|u|}^{2T-u}F^{11}(u,v)du$
$\displaystyle=-F^{11}|_{v=2T-u}-F^{11}|_{v=|u|}sgn(u)+\int_{|u|}^{2T-u}\partial_{u}F^{11}dv$
$\displaystyle=-F^{11}|_{v=2T-u}-F^{11}|_{v=|u|}sgn(u)-\int_{|u|}^{2T-u}\partial_{v}F^{12}dv+\int_{|u|}^{2T-u}|G_{1}|dv$
$\displaystyle=-F^{11}|_{v=2T-u}-F^{12}|_{v=2T-u}-F^{11}|_{v=|u|}sgn(u)+F^{12}|_{v=|u|}+\int_{|u|}^{2T-u}|G_{1}|dv$
$\displaystyle\lesssim(F^{11}+F^{12})|_{t=0}+\int_{|u|}^{2T-u}|G_{1}|dv.$

Integration in $u$ yields the desired estimate.

######  Lemma 3.2.

Under the assumption of Lemma3.1, we have

$\displaystyle\int_{0}^{T}\int_{0}^{\infty}(F^{11}F^{22}+F^{12}F^{21})(t,r)drdt$
$\displaystyle\lesssim\left(\int_{0}^{\infty}(F^{11}+F^{12})(0,r)dr+\int_{0}^{T}\int_{0}^{\infty}|G_{1}(t,r)|drdt\right)$
$\displaystyle\cdot\left(\int_{0}^{\infty}(F^{21}+F^{22})(0,r)dr+\int_{0}^{\infty}(F^{21}+F^{22})(T,r)dr+\int_{0}^{T}\int_{0}^{\infty}|G_{2}(t,r)|drdt\right).$

Proof Similary to the proof of Lemma3.1, for $\bar{v}\leq 2T-u$, we get

$\displaystyle\partial_{u}\int_{|u|}^{\bar{v}}F^{11}(u,v)dv+F^{12}(u,\bar{v})$
$\displaystyle\lesssim(F^{11}+F^{12})|_{u+\bar{v}=0}+\int_{|u|}^{\bar{v}}G_{1}dv.$

Thus

$\displaystyle
F^{21}(u,v)\partial_{u}\int_{|u|}^{\bar{v}}F^{11}(u,v)dv+F^{21}(u,\bar{v})F^{12}(u,\bar{v})$
$\displaystyle\lesssim
F^{21}(u,\bar{v})(F^{11}+F^{12})|_{u+\bar{v}=0}+\left(\int_{|u|}^{\bar{v}}G_{1}(u,v)dv\right)F^{21}(u,\bar{v}).$

By the second equation in (3.1), the first term is

$\displaystyle\partial_{u}\left[F^{21}(u,\bar{v})\int_{|u|}^{\bar{v}}F^{11}(u,v)dv\right]-\partial_{u}F^{21}(u,\bar{v})\int_{|u|}^{\bar{v}}F^{11}(u,v)dv$
$\displaystyle=\partial_{u}\left[F^{21}(u,\bar{v})\int_{|u|}^{\bar{v}}F^{11}(u,v)dv\right]-\partial_{\bar{v}}F^{22}(u,\bar{v})\int_{|u|}^{\bar{v}}F^{11}(u,v)dv$
$\displaystyle-G_{2}(u,\bar{v})\int_{|u|}^{\bar{v}}F^{11}(u,v)dv$
$\displaystyle=\partial_{u}\left[F^{21}(u,\bar{v})\int_{|u|}^{\bar{v}}F^{11}(u,v)dv\right]-\partial_{\bar{v}}\left[F^{22}(u,\bar{v})\int_{|u|}^{\bar{v}}F^{11}(u,v)dv\right]$
$\displaystyle+F^{22}(u,\bar{v})F^{11}(u,\bar{v})-G_{2}(u,\bar{v})\int_{|u|}^{\bar{v}}F^{11}(u,v)dv.$

Therefore, we get

$\displaystyle
F^{22}(u,\bar{v})F^{11}(u,\bar{v})+F^{21}(u,\bar{v})F^{12}(u,\bar{v})$
$\displaystyle\lesssim\partial_{\bar{v}}\left[F^{22}(u,\bar{v})\int_{|u|}^{\bar{v}}F^{11}(u,v)dv\right]-\partial_{u}\left[F^{21}(u,\bar{v})\int_{|u|}^{\bar{v}}F^{11}(u,v)dv\right]$
$\displaystyle+G_{2}(u,\bar{v})\int_{|u|}^{\bar{v}}F^{11}(u,v)dv+F^{21}(u,\bar{v})(F^{11}+F^{12})|_{u+\bar{v}=0}+\left(\int_{|u|}^{\bar{v}}G_{1}(u,v)dv\right)F^{21}(u,\bar{v}).$

Integration in the region $\left\\{0\leq u+\bar{v}\leq 2T,\ \bar{v}-u\geq
0\right\\}$ and using Lemma3.1 for the boundary estimate, we get the desired
conclusions.

##  4 Small energy implies regularity

In this section, we assume

$\displaystyle E_{0}\leq\epsilon_{1}$ (4.1)

is sufficiently small. Under this assumption, we shall prove the solution is
smooth. For that purpose, we only need to give the $H^{2}$ estimate of the
solution, see [4].

Let $\Psi=\Phi_{t}$. Differentiating the equation(2.1), we get

$\displaystyle\Psi_{tt}-\Psi_{rr}-\frac{\Psi_{r}}{r}$
$\displaystyle=4B(\Phi)(\Psi_{u},\Phi_{v})+4B(\Phi)(\Phi_{u},\Psi_{v})+4\Psi\cdot
B^{\prime}(\Phi)(\Phi_{u},\Phi_{v}),$ (4.2)

thus

$\displaystyle
r\Psi_{t}\cdot\left(\Psi_{tt}-\Psi_{rr}-\frac{\Psi_{r}}{r}\right)$
$\displaystyle=4r\Psi_{t}\cdot B(\Phi)(\Psi_{u},\Phi_{v})+4r\Psi_{t}\cdot
B(\Phi)(\Phi_{u},\Psi_{v})+4\Psi_{t}\cdot\left[\Psi\cdot
A^{\prime}(\Phi)(\Phi_{u},\Phi_{v})\right],$ $\displaystyle\Psi\cdot
B(\Phi)=0.$

Differentiating with respect to $t$ yields

$\displaystyle\Psi_{t}\cdot B(\Phi)=-\Psi\cdot\left[\Psi\cdot
A^{\prime}(\Phi)\right].$

So we get

$\displaystyle\partial_{t}\left[\frac{1}{2}(\Psi_{t}^{2}+\Psi_{r}^{2})r\right]-\partial_{r}\left[r\Psi_{t}\cdot\Psi_{r}\right]$
$\displaystyle=4r\Psi\cdot\left[\Psi\cdot
B^{\prime}(\Phi)\right](\Psi_{u},\Phi_{v})+4r\Psi\cdot\left[\Psi\cdot
B^{\prime}(\Phi)\right](\Phi_{u},\Psi_{v})+4r\Psi_{t}\cdot\left[\Psi\cdot
B^{\prime}(\Phi)\right](\Phi_{u},\Phi_{v})$ $\displaystyle=G_{1}.$ (4.3)

In a similar way, let $\hat{\Psi}=\Phi_{r}(t,r)$. We get

$\displaystyle\partial_{t}\left[\frac{1}{2}(\hat{\Psi}_{t}^{2}+\hat{\Psi}_{r}^{2}+\frac{\hat{\Psi}^{2}}{r^{2}})r\right]-\partial_{r}\left[r\hat{\Psi}_{t}\cdot\hat{\Psi}_{r}\right]$
$\displaystyle=4r\hat{\Psi}\cdot\left[\hat{\Psi}\cdot
B^{\prime}(\Phi)\right](\hat{\Psi}_{u},\Phi_{v})+4r\hat{\Psi}\cdot\left[\hat{\Psi}\cdot
B^{\prime}(\Phi)\right](\Phi_{u},\hat{\Psi}_{v})+4r\hat{\Psi}_{t}\cdot\left[\hat{\Psi}\cdot
B^{\prime}(\Phi)\right](\Phi_{u},\Phi_{v})$ $\displaystyle=\tilde{G}_{1}.$
(4.4)

Integrating in $t,r$, we get

$\displaystyle\sup_{0\leq t\leq
T}|D^{2}\Phi|^{2}_{L^{2}(\mathbb{R}^{2})}(t)\lesssim|\nabla^{2}u_{0}|^{2}_{L^{2}(\mathbb{R}^{2})}+|\nabla
u_{1}|^{2}_{L^{2}(\mathbb{R}^{2})}+\int_{0}^{T}\int_{0}^{\infty}(|\tilde{G}_{1}|+|G_{1}|)(t,r)drdt.$
(4.5)

The main purpose of this section is to prove

$\displaystyle\int_{0}^{T}\int_{0}^{\infty}(|\tilde{G}_{1}|+|G_{1}|)(t,r)drdt\lesssim\epsilon_{1}E_{1},$
(4.6)

where

$\displaystyle E_{1}=|\nabla^{2}u_{0}|^{2}_{L^{2}(\mathbb{R}^{2})}+|\nabla
u_{1}|^{2}_{L^{2}(\mathbb{R}^{2})}.$ (4.7)

For that purpose, we use an induction argument and first we assume

$\displaystyle\int_{0}^{T}\int_{0}^{\infty}(|\tilde{G}_{1}|+|G_{1}|)(t,r)drdt\lesssim
E_{1}.$ (4.8)

Then we get

$\displaystyle|D^{2}\Phi|^{2}_{L^{2}(\mathbb{R}^{2})}(t)\lesssim E_{1}.$ (4.9)

To obtain (4.6), we only estimate $G_{1}$. The estimate of $\tilde{G}_{1}$ can
be done in a similar way. By the expression of $G_{1}$, we get

$\displaystyle|G_{1}|$
$\displaystyle\lesssim|\Psi|^{2}(|\Psi_{u}||\Phi_{v}|+|\Phi_{u}||\Psi_{v}|)r+|\Psi_{t}|(|\Psi||\Phi_{u}||\Phi_{v}|)r$
$\displaystyle=g_{1}+g_{2}.$ (4.10)

We first estimate $g_{1}$. Let $\frac{1}{10}\leq\beta\leq\frac{1}{4}$, we have

$\displaystyle\int_{0}^{T}\int_{0}^{\infty}|g_{1}|drdt\lesssim\left[\int_{0}^{T}\int_{0}^{\infty}\left(|\Psi_{u}|^{2}|\Phi_{v}|^{2}+|\Phi_{u}|^{2}|\Psi_{v}|^{2}\right)r^{2-\beta}drdt\right]^{\frac{1}{2}}\left[\int_{0}^{T}\int_{0}^{\infty}|\Psi|^{4}r^{\beta}drdt\right]^{\frac{1}{2}}.$
(4.11)

By Sobolev-Hardy inequality,

$\displaystyle\int_{0}^{\infty}|\Psi|^{4}r^{\beta}dr\lesssim\left(\int_{0}^{\infty}|\Psi|^{2}r^{-\beta}dr\right)\left(\int_{0}^{\infty}|\Psi|^{2}rdr\right)^{\beta}\left(\int_{0}^{\infty}|\Psi_{r}|^{2}rdr\right)^{1-\beta},$
(4.12)

thus

$\displaystyle\int_{0}^{T}\int_{0}^{\infty}|g_{1}|drdt$
$\displaystyle\lesssim\left[\int_{0}^{T}\int_{0}^{\infty}\left(|\Psi_{u}|^{2}|\Phi_{v}|^{2}+|\Phi_{u}|^{2}|\Psi_{v}|^{2}\right)r^{2-\beta}drdt\right]^{\frac{1}{2}}\left[\int_{0}^{T}\int_{0}^{\infty}|\Psi|^{2}r^{-\beta}drdt\right]^{\frac{1}{2}}$

Take $\alpha=\beta$ in (2) and integrate for $0\leq r\leq+\infty,\ 0\leq t\leq
T$. Noting (2.19), we get

$\displaystyle\frac{1}{2}\int_{0}^{T}\int_{0}^{\infty}r^{-\beta}\left(|\Psi|^{2}+|\hat{\Psi}|^{2}\right)drdt=\left(\int_{0}^{\infty}r^{1-\beta}\Psi\hat{\Psi}\right)|_{0}^{T}\left(q^{2}_{0}+q^{2}_{1}\right)drdt$
$\displaystyle-\frac{1}{2}\int_{0}^{\infty}r^{-\beta+\sigma}Q_{0}q_{0}dr|_{0}^{T}$
$\displaystyle+\left[\frac{1}{2}\int_{0}^{\infty}q_{0}r^{\sigma-\alpha}\int_{0}^{r}\xi^{\alpha}A_{0}q_{1}d\xi-\sigma
q_{0}r^{\sigma-\alpha}\int_{0}^{r}\xi^{-\sigma-1}q_{0}(\xi)d\xi+\frac{1}{2}r^{\sigma-\alpha}A_{0}q_{0}Q_{0}drdt\right]$

By Hardy’s inequality and noting (2.10), the third term in the right hand of
side of the above equality can be absorbed by the left hand side, and by
Sobolev-hardy inequality the first and second term is bounded by

$\displaystyle\epsilon_{1}^{1-\beta/2}E^{\beta/2}_{1},$

thus

$\displaystyle\left(\int_{0}^{T}\int_{0}^{\infty}r^{-\beta}\left(|\Psi|^{2}+|\hat{\Psi}|^{2}\right)drdt\right)^{\frac{1}{2}}\lesssim\epsilon_{1}^{\frac{1}{2}-\frac{\beta}{4}}E_{1}^{\frac{\beta}{4}}.$
(4.14)

Rewrite (4.3) as

$\displaystyle\partial_{u}\left[r(\partial_{v}\Psi)^{2}\right]+\partial_{v}\left[r(\partial_{u}\Psi)^{2}\right]=2G_{1}.$
(4.15)

Take $\alpha=\beta$ in (2) and rewrite it as

$\displaystyle-\partial_{v}\left(r^{1-\beta}q_{0}q_{1}\right)+\partial_{v}\left(\frac{1}{2}r^{1-\beta}(q_{0}^{2}+q_{1}^{2})\right)-\frac{1}{2}\partial_{v}(r^{-\beta+\sigma}Q_{0}q_{0})$
$\displaystyle+\frac{1}{2}\partial_{v}(q_{1}r^{-\beta+\sigma}Q_{0})+\frac{\beta-\sigma+1}{4}\partial_{v}(r^{2\sigma-1-\beta}Q_{0}^{2})$
$\displaystyle-\partial_{u}\left(r^{1-\beta}q_{0}q_{1}\right)-\partial_{u}\left(\frac{1}{2}r^{1-\beta}(q_{0}^{2}+q_{1}^{2})\right)-\frac{1}{2}\partial_{u}(r^{-\beta+\sigma}Q_{0}q_{0})$
$\displaystyle-\frac{1}{2}\partial_{u}(q_{1}r^{-\beta+\sigma}Q_{0})-\frac{\beta-\sigma+1}{4}\partial_{u}(r^{2\sigma-1-\beta}Q_{0}^{2})$
$\displaystyle=\frac{q_{0}}{2r^{\beta-\sigma}}\int_{0}^{r}\xi^{-\sigma}A_{0}q_{1}d\xi-\frac{\sigma
q_{0}}{2r^{\beta-\sigma}}\int_{0}^{r}\xi^{-\sigma-1}q_{0}(\xi)d\xi+\frac{1}{2}r^{-\beta+\sigma}A_{0}q_{0}Q_{0}$
$\displaystyle-\frac{(1+\beta-2\sigma)(\beta-\sigma+1)}{4}r^{2\sigma-2-\alpha}Q_{0}^{2}-\frac{\beta}{2}\frac{q_{0}^{2}+q_{1}^{2}}{r^{\alpha}}$
. (4.16)

$\displaystyle\partial_{v}\left(r^{1-\beta}\frac{1}{2}\left(q_{0}-q_{1}\right)^{2}+\frac{1}{2}r^{-\beta+\sigma}Q_{0}\left(q_{1}-q_{0}\right)+\frac{\beta-\sigma+1}{4}r^{2\sigma-1-\beta}Q_{0}^{2}\right)$
$\displaystyle-\partial_{u}\left(r^{1-\beta}\frac{1}{2}\left(q_{0}+q_{1}\right)^{2}+\frac{1}{2}r^{-\beta+\sigma}Q_{0}\left(q_{1}+q_{0}\right)+\frac{\beta-\sigma+1}{4}r^{2\sigma-1-\beta}Q_{0}^{2}\right)$
$\displaystyle=G_{\beta},$ (4.17)

where

$\displaystyle G_{\beta}:=$
$\displaystyle\frac{q_{0}}{2r^{\beta-\sigma}}\int_{0}^{r}\xi^{-\sigma}A_{0}q_{1}d\xi-\frac{\sigma
q_{0}}{2r^{\beta-\sigma}}\int_{0}^{r}\xi^{-\sigma-1}q_{0}(\xi)d\xi+\frac{1}{2}r^{-\beta+\sigma}A_{0}q_{0}Q_{0}$
$\displaystyle-\frac{(1+\beta-2\sigma)(\beta-\sigma+1)}{4}r^{2\sigma-2-\alpha}Q_{0}^{2}-\frac{\beta}{2}\frac{q_{0}^{2}+q_{1}^{2}}{r^{\alpha}}.$
(4.18)

We integrate above balace law over space time, and take the $\sigma$
sufficiently small. By hardy inequality, we are arrive at the following
estimates

$\displaystyle\frac{(1+\beta-2\sigma)(\beta-\sigma+1)}{4}r^{2\sigma-2-\alpha}Q_{0}^{2}+\frac{\beta}{2}\frac{q_{0}^{2}+q_{1}^{2}}{r^{\alpha}}$
$\displaystyle\leq
2\left(\frac{(1+\beta-2\sigma)(\beta-\sigma+1)}{4}r^{2\sigma-2-\alpha}Q_{0}^{2}+\frac{\beta}{2}\frac{q_{0}^{2}+q_{1}^{2}}{r^{\alpha}}+\frac{\sigma
q_{0}}{2r^{\beta-\sigma}}\int_{0}^{r}\xi^{-\sigma-1}q_{0}(\xi)d\xi\right)$
$\displaystyle\leq 2\left(\frac{\sigma
q_{0}}{2r^{\beta-\sigma}}\int_{0}^{r}\xi^{-\sigma}A_{0}q_{0}(\xi)d\xi+\frac{1}{2}r^{-\beta+\sigma}A_{0}q_{0}Q_{0}\right)+C\epsilon^{1-\frac{\beta}{2}}_{1}E^{\frac{\beta}{2}}_{1}$

Noting that (2.10) and energy small, the first two terms on the right can be
absorbed by the left hand side, therefore

$\displaystyle\int_{0}^{T}\int_{0}^{\infty}|G_{\beta}|(t,r)drdt\lesssim\epsilon_{1}^{1-\frac{\beta}{2}}E_{1}^{\frac{\beta}{2}}.$
(4.20)

Noting that

$\displaystyle
r^{1-\beta}\lvert\Phi_{u}\rvert^{2}=\frac{1}{4}r^{1-\beta}\left(q_{0}-q_{1}\right)^{2}$
$\displaystyle\lesssim
r^{1-\beta}\frac{1}{2}\left(q_{0}-q_{1}\right)^{2}+\frac{1}{2}r^{-\beta+\sigma}Q_{0}\left(q_{1}-q_{0}\right)+\frac{\beta-\sigma+1}{4}r^{2\sigma-1-\beta}Q_{0}^{2}$
(4.21)

and

$\displaystyle
r^{1-\beta}\lvert\Phi_{v}\rvert^{2}=\frac{1}{4}r^{1-\beta}\left(q_{0}-q_{1}\right)^{2}$
$\displaystyle\lesssim
r^{1-\beta}\frac{1}{2}\left(q_{0}+q_{1}\right)^{2}+\frac{1}{2}r^{-\beta+\sigma}Q_{0}\left(q_{1}+q_{0}\right)+\frac{\beta-\sigma+1}{4}r^{2\sigma-1-\beta}Q_{0}^{2}$
(4.22)

We apply Lemma 3.2 to get

$\displaystyle\int_{0}^{T}\int_{0}^{\infty}r^{2-\beta}\left[(\partial_{u}\Phi)^{2}(\partial_{v}\Psi)^{2}+(\partial_{u}\Psi)^{2}(\partial_{v}\Phi)^{2}\right]drdt\lesssim\epsilon_{1}^{1-\frac{\beta}{2}}E_{1}^{1+\frac{\beta}{2}}.$
(4.23)

So we get

$\displaystyle\int_{0}^{T}\int_{0}^{\infty}|g_{1}|drdt\lesssim\epsilon_{1}E_{1}.$
(4.24)

Now we estimate $g_{2}$.

$\displaystyle\int_{0}^{T}\int_{0}^{\infty}|g_{2}|drdt\lesssim\int_{0}^{T}\left(\int_{0}^{\infty}|\Psi_{t}|^{2}rdr\right)^{\frac{1}{2}}\left(\int_{0}^{\infty}|\Psi|^{2}|\Phi_{u}|^{2}|\Phi_{v}|^{2}rdr\right)^{\frac{1}{2}}dt$
$\displaystyle\lesssim
E_{1}^{\frac{1}{2}}\int_{0}^{T}\left(\int_{0}^{\infty}|\Psi|^{2}|\Phi_{u}|^{2}|\Phi_{v}|^{2}rdr\right)^{\frac{1}{2}}dt.$
(4.25)

We have

$\displaystyle\int_{0}^{\infty}|\Psi|^{2}|\Phi_{u}|^{2}|\Phi_{v}|^{2}rdr=\int_{0}^{\infty}|\Phi_{u}|^{2}|\Phi_{v}|^{2}d\int_{0}^{r}|\Psi|^{2}\lambda
d\lambda$ $\displaystyle=-2\int_{0}^{\infty}\int_{0}^{r}|\Psi|^{2}\lambda
d\lambda\left[(\Phi_{u}\cdot\hat{\Psi}_{u})|\Phi_{v}|^{2}+|\Phi_{u}|^{2}(\Phi_{v}\cdot\hat{\Psi}_{v})\right]$
$\displaystyle\leq
2\int_{0}^{\infty}dr\int_{0}^{r}|\Psi|^{2}\lambda^{-\beta}d\lambda\left[r^{1+\beta}|\Phi_{u}||\Phi_{v}|(|\hat{\Psi}_{u}||\Phi_{v}|+|\hat{\Psi}_{v}||\Phi_{u}|)\right]$
$\displaystyle\leq
2\left(\int_{0}^{\infty}|\Psi|^{2}r^{-\beta}dr\right)\left(\int_{0}^{\infty}(|\hat{\Psi}_{u}|^{2}|\Phi_{v}|^{2}+|\hat{\Psi}_{v}|^{2}|\Phi_{u}|^{2})r^{2-\beta}dr\right)^{\frac{1}{2}}$
$\displaystyle\left(\int_{0}^{\infty}|\Phi_{u}|^{2}|\Phi_{v}|^{2}r^{3\beta}dr\right)^{\frac{1}{2}}.$

Therefore we get

$\displaystyle\int_{0}^{T}\int_{0}^{\infty}|g_{2}|drdt\lesssim
E_{1}^{\frac{1}{2}}\left(\int_{0}^{T}\int_{0}^{\infty}|\Psi|^{2}r^{-\beta}drdt\right)^{\frac{1}{2}}$
$\displaystyle\left(\int_{0}^{T}\int_{0}^{\infty}r^{2-\beta}(|\hat{\Psi}_{u}|^{2}|\Phi_{v}|^{2}+|\hat{\Psi}_{v}|^{2}|\Phi_{u}|^{2})drdt\right)^{\frac{1}{4}}$
$\displaystyle\left(\int_{0}^{T}\int_{0}^{\infty}|\Phi_{u}|^{2}|\Phi_{v}|^{2}r^{3\beta}drdt\right)^{\frac{1}{4}}.$
(4.26)

The first term in the right hand side can be estimated by (4.14), and the
second term can be estimated in a way similar to (4.23). Thus, it remains to
prove

$\displaystyle\int_{0}^{T}\int_{0}^{\infty}|\Phi_{u}|^{2}|\Phi_{v}|^{2}r^{3\beta}drdt\lesssim\epsilon_{1}^{1+\frac{3\beta}{2}}E_{1}^{1-\frac{3\beta}{2}}.$

We have by (4.17) with $\beta$ replaced by $\alpha$

$\displaystyle\partial_{v}\left(r^{1-\alpha}\frac{1}{2}\left(q_{0}-q_{1}\right)^{2}+\frac{1}{2}r^{-\alpha+\sigma}Q_{0}\left(q_{1}-q_{0}\right)+\frac{\alpha-\sigma+1}{4}r^{2\sigma-1-\alpha}Q_{0}^{2}\right)$
$\displaystyle-\partial_{u}\left(r^{1-\alpha}\frac{1}{2}\left(q_{0}+q_{1}\right)^{2}+\frac{1}{2}r^{-\alpha+\sigma}Q_{0}\left(q_{1}+q_{0}\right)+\frac{\alpha-\sigma+1}{4}r^{2\sigma-1-\alpha}Q_{0}^{2}\right)$
$\displaystyle=G_{\alpha},$ (4.27)

We have

$\displaystyle\Phi_{t}\cdot\left[\Phi_{tt}-\Phi_{rr}-\frac{\Phi_{r}}{r}\right]=0,$

so we get

$\displaystyle\partial_{u}\left[r^{1-\alpha}(\partial_{v}\Phi)^{2}\right]+\partial_{v}\left[r^{1-\alpha}(\partial_{u}\Phi)^{2}\right]=\hat{G}_{\alpha},$

where

$\displaystyle\hat{G}_{\alpha}=(1-\alpha)r^{-\alpha}\left[(\partial_{v}\Phi)^{2}-(\partial_{u}\Phi)^{2}\right].$

Thus, we have

$\displaystyle\int_{0}^{T}\int_{0}^{\infty}\left(|\hat{G}_{\alpha}|+|G_{\alpha}|\right)drdt\lesssim\epsilon_{1}^{1-\frac{\alpha}{2}}E_{1}^{\frac{\alpha}{2}}.$

So we apply Lemma3.2 to get

$\displaystyle\int_{0}^{T}\int_{0}^{\infty}r^{2-2\alpha}(\partial_{v}\Phi)^{2}(\partial_{u}\Phi)^{2}drdt\lesssim\epsilon_{1}^{2-\alpha}E_{1}^{\alpha}.$

Take $2-2\alpha=3\beta$, we get

$\displaystyle\int_{0}^{T}\int_{0}^{\infty}r^{3\beta}(\partial_{v}\Phi)^{2}(\partial_{u}\Phi)^{2}drdt\lesssim\epsilon_{1}^{1+\frac{3\beta}{2}}E_{1}^{1-\frac{3\beta}{2}}.$

This completes the proof of the regularity with small energy.

######  Remark 4.1.

This paper was published in Chinese Annals of Mathematics in 2022,
unfortunately there is a miss calculation in the main equality. The present
version correct this mistakes.

Acknowledgement The author is supported by Key Laboratory of Mathematics for
Nonlinear Sciences (Fudan University), Ministry of Education of China,
P.R.China. Shanghai Key Laboratory for Contemporary Applied Mathematics,
School of Mathematical Sciences, Fudan University, P.R. China, and by Shanghai
Science and Technology Program [Project No. 21JC1400600].

## References

  * [1] Christodoulou, D. and Tahvildar-Zadeh, A.S., On the regularity of spherically symmetric wave maps, Comm. Pure Appl. Math, 46, 1993, 1041–1091.
  * [2] Geba, D.A. and Grillakis M.G., An introduction to the theory of wave maps and related geometric problems, World Scientific Publishing Co., 2016.
  * [3] Gu, Chao-Hao, On the Cauchy problem for harmonic maps defined on two-dimensional Minkovski space, Commun. Pure Appl. Math. ,33, (1980), 727-737.
  * [4] Shatah, J. and Struwe, M., Geometric wave equations, Courant Lecture Notes in Mathematics, 2, 1998.
  * [5] Shatah, J. and Struwe, M., The Cauchy problem for wave maps, Int. Math. Res. Not., 11, 2002, 555–571.
  * [6] Struwe, M., Radially symmetric wave maps from (1 + 2)-dimensional Minkowski space to the sphere, Math. Z., 242, 2002, 407–414.
  * [7] Struwe, M., Radially symmetric wave maps from (1 + 2)-dimensional Minkowski space to general targets, Cal. Var., 16, 2003, 431–437.
  * [8] Tao, T., Nonlinear dispersive equations, CBMS Regional Conference Series in Math., 106, 2006.


rather they are high-order tensor contractions. Furthermore, because of energy
and momentum conservation in these nonlinear interactions, it will rarely be
possible to realize fully arbitrary high-order tensor contractions or
convolutions. While PAT and PNNs can facilitate training physical computations
based on these exotic operations (this was one scientific motivation for using
nonlinear optics as a first demonstration of PNNs here), it is not yet clear
how those PNNs should be designed to be useful, since machine learning based
on their underlying operations is not well-studied. In summary, while single-
mode, nanophotonic nonlinear optics presents many promising features for PNNs,
we believe that significant effort will be required to understand how to
utilize the exotic nonlinear optical operations in practical tasks. In
contrast, other systems (such as multimode waves, coupled oscillators or
electrical networks) discussed subsequently have a more direct connection to
standard machine learning operations, so we expect it will be easier in the
short term to apply these systems as PNNs, using standard machine learning
techniques, to practical problems.

###  3.3 Self-simulation analysis of multimode wave propagation

Multimode wave propagation is an appealing platform for physical neural
networks for similar reasons that motivate its use in physical reservoir
computing [27, 28]. In the linear regime, multimode wave propagation permits
very high-dimensional random matrix-vector multiplications with potentially
very low-energy [27]. If the medium light propagates through can be
controlled, then these matrix-vector operations will likewise be controllable,
making for a very direct connection to standard deep neural network
operations. In the nonlinear regime, the complexity of pulse propagation can
grow significantly, allowing for PNNs with much greater self-simulation
advantages. This logic has motivated using multimode nonlinear wave
propagation for reservoir computing [28].

####  3.3.1 Linear propagation

All linear multimode wave propagation can be modelled as a matrix-vector
operation, where a transmission matrix $T$ operates on an input vector of
complex mode amplitudes $\vec{x}$ [29]. Thus, the number of digital operations
required to simulate a linear multimode wave system scales as $MM^{\prime}$,
where $M$ is the number input modes and $M^{\prime}$ the number of measured
output modes. Since the energy required to implement the wave propagation
scales linearly with $M$ (assuming $M^{\prime}=M$), and the energy required to
perform the matrix-vector multiplication scales as $MR$, where $R$ is the rank
of the transmission matrix, the self-simulation advantage scales linearly with
the number of modes.

As an example, multimode fibers can support $\sim 5\times 10^{6}$ modes, and
fibers of $\sim 10$ cm length can exhibit all-to-all random coupling [30]. An
experiment involving a 1 mW beam, a 30 kHz digital micromirror device with 5
million pixels, and a high-speed 1 Mpixel camera would consume approximately
10 mW + 1 W + 5 W $\approx$ 6 W and performs, in effect, $(5\times
10^{6})^{2}\times 30\times 10^{3}\text{Hz}\sim 8\times 10^{17}$ operations per
second, corresponding to about 10 aJ/operation. In comparison, the GTX GPU
(DGX Superpod) requires $10^{6}$ ($5\times 10^{6}$) times more energy, and
$10^{4}$ (300) times more time. Even if such an operation has a limited number
of trainable parameters, large random matrix operations appear to be useful
for a variety of applications [31, 27]. As components within PNNs, they could
be integrated within conventional neural architectures. Their parameters could
be trained using gradient descent alongside others using PAT. Even in
implementations without trainable parameters, PAT would also allow
backpropagation through their operations, to train preceding or subsequent
trainable operations as part of hybrid digital-physical architectures.
Nonetheless, it is obvious that having controllable physical parameters to
permit arbitrary, or nearly-arbitrary matrix-vector operations would give rise
to a stronger connection with the format of standard deep neural networks, and
therefore probably to a faster route to good performance on practical tasks.

####  3.3.2 Nonlinear propagation

In the nonlinear regime, multimode wave propagation is considerably more
complex, and the self-simulation advantages can be very high. While such large
self-simulation advantages recommend complex nonlinear wave systems (and many-
body nonlinear dynamical systems in general), we caution that these large
numbers should be taken with significant skepticism with respect to their
translation to general-purpose advantages across a wide range of practical
problems. While the performance of multimode nonlinear wave propagation-based
reservoir computers suggests that some practical advantages may be obtained
[28], it is still challenging to finely control multimode nonlinear waves and
the main operations that contribute to the large self-simulation advantages
are exotic to modern neural networks. Overall, much more work is necessary.

Additionally, systems with extremely large self-simulations advantages should
be evaluated with caution because it may be possible to dramatically
accelerate the simulations of these systems not only by more efficient
simulation algorithms, but also by data-driven simulations, such as in Ref.
[32]. This will be especially likely if the simulations need only to consider
a relatively constrained set of possible inputs, rather than a fully general
simulation, or when interactions are subject to symmetry constraints.

With these caveats noted, we now consider multimode nonlinear optical wave
propagation. The standard model for multimode Kerr media scales with the
number of guided modes $M$ as $M^{4}$ [33]. Although this scaling can be
significantly alleviated if symmetry assumptions can be applied [34], in
fibers with strong disorder (which implements all-to-all linear coupling as
considered above) these assumptions will be broken and the $M^{4}$ scaling
will hold. In highly-multimode fibers, more efficient simulations may be
performed by simulating light propagation using a 3+1-dimensional split-step
Fourier algorithm, similar to the one considered in the previous section [35].
However, these simulations still scale at least as fast as
$M^{2}N_{\omega}/E_{p}L_{b}$ in the general case, where $M$ is the number of
guided modes, $E_{p}$ is the pulse energy, $N_{\omega}$ is the number of
frequency points, and $L_{b}$ is the beat length for interference between
involved modes. Since $L_{b}\sim 1/\Delta_{n}$ grows with the number of modes
in the fiber, the effective scaling is roughly $M^{\alpha}N_{\omega}/E_{p}$,
where $\alpha$ is a scaling parameter at least larger than $2$.

In the former experiment, assuming a 1 W pulsed laser source with 31 kHz
repetition rate, $32\text{\,}\mathrm{\SIUnitSymbolMicro J}$ pulses would
produce complex supercontinuum in even highly-multimode fibers such as the
$5\times 10^{6}$ one considered above. To simulate such pulse propagation,
$N_{x}>3000$ spatial grid points ($\approx\sqrt{5\times 10^{6}}$, the best-
case number for a fiber with $5\times 10^{6}$ modes) and $N_{\omega}>10^{4}$
($\approx$ 100 THz $\times$ 100 ps, the expected spectral bandwidth times a
minimal estimate of the temporal window required). Assuming an core-cladding
index difference of 0.01 and a 532 nm center wavelength, the beat length of
$\Delta_{n}/\lambda$ is $53\text{\,}\mathrm{\SIUnitSymbolMicro m}$, suggesting
a best-case step-size of $\sim$ $25\text{\,}\mathrm{\SIUnitSymbolMicro m}$.
Then, the number of operations required to simulate a 100-cm fiber evolution
would be approximately:

$[N_{x}\log_{2}(N_{x})]^{2}N_{\omega}\log_{2}N_{\omega}L/L_{\text{step}}=6.4\times
10^{18}\text{ operations.}$ (37)

Since these operations are, in effect, physically executed at a 31 kHz rate
(limited by the DMD update), the effective operation rate is $2\times 10^{23}$
operations/s. Assuming the same components as described above and a 1% laser
wall-plug efficiency, this corresponds to $106\text{\,}\mathrm{W}$/$2\times
10^{23}=5\times 10^{-22}$ J/operation. This number is astoundingly small: the
physics is $10^{10}$ ($10^{11}$) more energy efficient than a digital
simulation using the RTX GPU (DGX Superpod), and executes 105 (2000) times
faster.

###  3.4 Multimode nonlinear oscillators/oscillator networks

Networks of coupled oscillators can be simulated using a set of differential
equations like 29, although for bulk multimode oscillators in which the modes
overlap significantly in space it will usually be more efficient to instead
simulate partial differential equations for the time-evolution of the bulk
field in $(x,y,z)$. As an example of a locally-connected network of localized
nonlinear oscillators, we consider nanoscale spintronic oscillators, such as
those used for reservoir computing in Ref. [10]. In that work, they claim that
the power consumption per oscillator could be as low as
$1\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{W}$, with resonance frequencies
into the 10s of GHz [36].

If we assume that nonlinearity in the oscillators is strictly local, the
number of operations required to simulate a network of coupled spintronic
oscillators is then:

$N_{\text{ops}}\approx\Delta T\Delta f[NM+N]$ (38)

where $\Delta T$ is the simulated time, $\Delta f$ is the relevant frequency
bandwidth, $N$ is the number of oscillators, $M$ is the number of nearest
neighbours each oscillator is (linearly) coupled to.

Given these assumptions, the coupled oscillator network follows an equation
that closely resembles a typical deep residual neural network [37]:

$\dot{x}_{i}(t)=\sum_{j=1}^{N}a_{ij}(t)x_{j}(t)+f_{i}(t)+\sigma(x_{i}(t)),$
(39)

where $\sigma$ is the local nonlinearity of each oscillator. As is noted in
neural ordinary differential equations [38], integrating the above equation
using, for example, an Euler method, leads to an update equation:

$x_{i}[t+1]=x_{i}[t]+\Delta
t(f_{i}[t]+\sum_{j=1}^{N}a_{ij}[t]x_{j}[t]+\sigma(x_{i}[t])),$ (40)

which is just the layer-by-layer propagation through a deep residual neural
network.

Let us assume a 10 GHz resonance, such that the maximum relevant frequency in
a simulation is 20 GHz. We’ll consider $N=100$ oscillators with a coupling
range of $M=100$ and assume that nonlinearity can be implemented by 1
operation per oscillator per step.

For a time of 1 ns, the number of operations is then:

$N_{\text{ops}}\approx\Delta T\Delta f[NM+N]=1\text{ ns}\times 20\text{
GHz}\times[100\times 100+100]\approx 2\times 10^{5}$ (41)

This costs $2\text{\,}\mathrm{\SIUnitSymbolMicro J}$
($7.9\text{\,}\mathrm{\SIUnitSymbolMicro J}$) and takes 6.3 ns (72 ps) on the
GTX GPU (DGX Superpod) respectively.

If we assume that each oscillator requires 2.5 times its nominal power
requirements to allow it to efficiently couple with 100 neighbours, each
oscillator requires $2.5\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{W}$ and
the total energy required physically is 250 fJ. However, the input and output
voltages signals will in general require conversion from/to digital format. As
discussed in the previous sections, a cost of 1 pJ/sample for each DAC and ADC
operation at 20 GHz is plausible [19]. For 100 inputs and 100 outputs with 20
GHz $\times$ 1 ns = 20 samples each, the DAC and ADC costs are each 100
$\times$ 20 $\times$ 1 pJ = 2 nJ. Thus, the total energy cost to evolve the
physical system is about 4 nJ.

Thus, the self simulation advantage of the network is about $10^{3}$ for
energy, and between about 0.1 to 10 for speed. While this seems relatively
small compared to previous subsections, it is straightforward to increase.

In this example, the energy required for the oscillators themselves was
negligible compared to the DAC and ADC costs. Thus, we can increase the
complexity of operations the device performs by increasing the ‘hidden’
oscillators. Let us now consider a network of $N=10^{5}$ coupled oscillators,
each coupled to $M=100$ neighbours, and of which only $N_{\text{in}}=100$ and
$N_{\text{out}}=100$ are used for input and output. The rest of the
oscillators we imagine to be used to cause and control the physical
transformation. Their injection currents might be trainable parameters which
can be changed at a much slower rate during training, and otherwise maintained
in-placed during inference (and therefore not incurring DAC/ADC costs each
inference). If we assume that
$10\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{W}$ is required to maintain
each of these “parameter” oscillators, then the total energy for oscillators
over 1 ns is about 1 nJ. This is still less than the energy required for DAC
and ADC.

For the larger oscillator network, the number of operations is however much
larger, $2\times 10^{8}$. Consequently, the calculation costs $10^{5}$ times
more than the first example: it costs 2 mJ (8 mJ) and takes
$6.3\text{\,}\mathrm{\SIUnitSymbolMicro s}$ (71 ns) on the GTX GPU (DGX
Superpod) respectively. And, as a result, the larger oscillator network’s SSAs
are around $10^{6}$ to $10^{7}$ for energy and 70 to 6000 for speed.

We note that greater speed advantages are possible if one increases $N$ while
holding $N\Delta T$ constant. When this is done, the SSA increases linearly
with $N$. Of course, there will be a minimum useful $\Delta T$, and as the
number of oscillators increases, it will usually be necessary to increase this
time so that distant ones can meaningfully interact. Thus, how much spatial
parallelism/scaling can be exploited is unclear and will depend on the system
and the required transformations.

Finally, we note that this section has neglected energy costs associated with
the controllable linear coupling. We have inherently assumed a technique that
requires no additional power to maintain, such as bi- or multi-stable circuit
elements. Even though these devices may require power to adjust, their
influence on the cost of performing inference could be, ideally, close to
zero. We will expand on this point at the end of the next sub-section, but
obviously if connections require active power to maintain, even very small
amounts, it could rapidly become non-negligible due to the enormous number of
connections involved.

###  3.5 Nonlinear electrical circuit dynamics/analog transistor networks

To describe the transient analog nonlinear dynamics of a network of
transistors, one may use a nonlinear circuit simulation strategy, solving a
system of nonlinear differential equations of at least size $N_{T}N_{d}$,
where $N_{T}$ is the number of transistors, and $N_{d}$ is the number of
dimensions ($\sim$nodes) in the equivalent circuit model of the transistor.
For accurate large-signal analog models such as the EKV model, $N_{d}\sim 10$.
To solve for the dynamics of a nonlinear network of transistors, a common
technique is to numerically integrate the differential equations that describe
the network by solving at each time-step the nonlinear system of equations
using a nonlinear solver such as Newton-Raphson. Alternatively, if we can
approximate all nonlinearities by low-order Taylor-series approximations, one
may instead integrate a system of coupled nonlinear differential equations
directly. Since the nonlinearities in electrical circuits are typically
strong, including strong saturation and exponential dependencies, the latter
approach is not often used.

The standard approach scales in complexity as $\mathcal{O}(n^{1.3-1.8}$) [39,
40, 41], with the most complex step due to the linear system solving required
for each iteration of the Newton-Raphson method. For simulating very large
circuits, there are methods to improve this scaling depending on the topology
and details of the circuit [40, 41]. However, these methods sacrifice
accuracy, and depend on the specific circuit topology.

Even though it is less accurate, integrating a system of equations similar to
29 allows more straightforward insight into the circuit simulation scaling and
the operations it effectively performs. Assuming integration of the nonlinear
system of differential equations, the number of operations to simulate the
transistor network is nominally:

$N_{\text{ops}}\approx\Delta T\Delta
f[N_{T}N_{d}M_{1}+N_{T}N_{d}^{2}M_{2}^{2}+...+N_{T}N_{d}^{n}M_{n}^{n}],$ (42)

where $\Delta T$ = time of simulation, $\Delta f$ = bandwidth of simulation
(the inverse of the time-step, roughly the maximum frequency expected to be
encountered in the dynamics), $N_{d}$ = number of degrees of freedom needed to
describe each transistor $\sim$ 10-20, $n$ = highest polynomial order of the
nonlinear coupling needed to describe the operating regime, and $M_{i}$ =
range of transistor-transistor coupling (i.e., the number of other transistors
each transistor is non-negligibly coupled to) for each order of nonlinearity.

Although it may not be true when nano-sized transistors are operated near one
another at radio frequencies, we will assume that the nonlinear terms are
local to each transistor, so $M_{i}=1$ for $i$ larger than 1. In this case,

$N_{\text{ops}}\approx\Delta T\Delta
f[N_{T}N_{d}M_{1}+N_{T}N_{d}C_{\text{nl}}],$ (43)

where we have assumed that a more efficient method is used to calculate the
nonlinearity at each transistor, such that the scaling is a constant factor
larger than the number of degrees of freedom required to simulate each
transistor. We have done this because there are obviously much more efficient
methods for the transistor nonlinearity than the full ODE expansion in Eqn.
29, and we would thus grossly overestimate the simulation cost. If we let
$C_{\text{nl}}=N_{d}$, it is in some sense equivalent to assuming that the
calculation of the nonlinearity is about as efficient as if the nonlinearity
were approximated by terms up to quadratic. Provided $C_{\text{nl}}$ is not
much larger than $M_{1}$, its choice is not significant in our final estimate,
since the nonlinear part of the simulation will be comparable to or smaller
than the linear coupling part.

Before considering a concrete example, it is worthwhile to emphasize that,
again, with these assumptions, and assuming that the internal complexities of
the transistor are abstracted away, the physical system follows an equation
for which each numerical step resembles propagation through a residual layer:

$x_{i}[t+1]=x_{i}[t]+\Delta
t(f_{i}[t]+\sum_{j=1}^{N}a_{ij}[t]x_{j}[t]+\sigma(x_{i}[t])).$ (44)

Although the transistor network differs in some details from the oscillator
networks just considered, they both ultimately resemble typical neural
networks provided their linear coupling is controllable, and their
nonlinearity is relatively local. Since simulating the equivalent circuit for
a transistor in the nonlinear, RF-regime is significantly harder than what we
assumed for the nonlinear oscillators, the transistor networks may have a
modest constant-factor advantage in self-simulation over nonlinear oscillators
(around 10-20 times). While PNNs could make use of such physics beyond Eqn.
44, it is difficult to estimate how much computational benefit such higher-
order physics would confer. The main takeaway is just that, up to some
approximations and abstractions, networks of coupled transistors (or diodes)
perform operations as they evolve that are very similar to those in
conventional neural networks. Therefore, like the oscillator networks in the
previous section, we feel relatively more confident in asserting that PNNs
based on these systems may ultimately be able to approach their self-
simulation advantage for tasks of practical interest (i.e., that computational
advantages for practical tasks could be as high as $\sim$ 1$\%$ of the self-
simulation advantages).

As an example, let us assume an electronic network of typical modern
transistors. Modern CPUs/GPUs require 10 W per billion transistors ($\sim$ 100
nW per transistor, on average) [42], but in part this power requirement relies
on efficient use of transistors such that all are not on all the time, so the
power consumption per transistor during its operation is much higher. Clock
rates of up to few GHz are typical. In typical digital usage, these
transistors are operated in constrained safe-operation regimes. However, as
analog devices, these same transistors are more complex, especially if driven
outside these typical regimes, for high-frequency analog signals beyond the
small-signal regime [43, 44]. We will therefore assume the transistors
dissipate roughly 1000 nW each, and can be operated at digital clock cycles of
a few GHz, such that their nonlinear analog dynamics can be as fast as 20 GHz
or so. In addition, we will typically want to use digital-analog conversion
for inputs. As is discussed in the last section, it is ideal if the majority
of the circuit’s parameters are stored in-place. We will therefore assume
there are $N_{T}$ simultaneous, time-dependent inputs to and outputs from the
circuit, requiring $N_{T}$ DAC and ADCs.

Let us assume $N_{T}$ = 100 transistors described by a $N_{d}=10$ equivalent
circuit, an evolution time of 10 ns, 5 GHz bandwidth inputs (so that a maximum
frequency of 20 GHz occurs in the simulation) and a coupling range of $M=100$
neighbouring transistors. This leads to a physical energy cost of:

$N_{T}\times P_{\text{{transistor}}}\times\Delta T=100\times 100\text{
nW}\times 10\text{ ns}=1\text{ pJ}$ (45)

for the transistors themselves, and

$N_{T}\times 2\times\Delta Tf_{\text{ s}}=100\times 10\text{ ns}\times 1\text{
pJ}\times 5\text{ GHz}=10\text{ nJ}$ (46)

for the DAC and ADC, where $f_{\text{s}}$ is the input and output sampling
frequency. Arguably, the transistor power requirement can be neglected, both
because it is much smaller than any DAC/ADC cost, but also because the signals
from the DAC provide the power required to drive the transistors. We will
however not do this, since we are interested in order-of-magnitude estimates
that are preferably conservative.

Meanwhile, simulations require:

$N_{\text{ops}}\approx\Delta T\Delta
f[N_{T}N_{d}M_{1}+N_{T}N_{d}C_{\text{nl}}]=10\text{ ns}\times 20\text{
GHz}\times[100\times 10\times 100+100\times 10\times 10]\approx 10^{7}$ (47)

As a result, the digital simulation, which consists effectively of locally-
structured matrix-vector multiplications and 3rd-order tensor contractions,
takes 0.1 mJ (0.4 mJ) and 310 ns (3.6 ns) on the GTX GPU (DGX Superpod).

From this, we see that a small ($N_{T}=100$) analog transistor network has a
self-simulation advantage of around $10^{4}$ to $10^{5}$ for energy, and
between 0.4 and 30 for speed. To increase this advantage (and thus,
speculatively, the computational capabilities of a PNN based on the analog
transistor network), it is clear that establishing longer-range coupling $M$
would be desirable.

Much like the previous section, we could increase the speed advantage simply
by increasing the number of transistors. For example, with $N_{T}=10^{3}$, and
$\Delta T=1$ ns, the number of operations required for simulation is the same
as the example above, and the total energy required to run the network for
$\Delta T$ is the same (since $N_{T}\Delta T$ was held constant). However, in
this case the electrical circuit has a speed advantage of 310 (4) over the GTX
GPU (DGX Superpod). Obviously there is a limit to how much this can be
exploited, since the evolution time of the circuit cannot be made arbitrarily
short if we require non-trivial physical transformations and long-range
coupling. Nonetheless, this illustrates how adding parallel physical processes
to a physical system can significantly enhance its potential computational
power.

Just like the previous section, we could also increase the self-simulation
advantages by adding additional “parameter” transistors that do not require
fast input and output, and thus do not add to the DAC and ADC costs, but
instead add only a smaller cost associated with providing enough total power
to compensate for the dissipation of the additional transistors. If
maintaining each additional such transistor requires an additional
$1\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{W}$, one could increase the
self-simulation advantages by nearly $10^{6}$ by adding $10^{6}$ such
transistors before the cost of these additional transistors was comparable to
the DAC/ADC cost.

Finally, we note that the above discussion has neglected the details of
realizing controllable linear coupling between transistors (or, in the earlier
section, oscillators). Our findings should hold for coupling techniques that
require no energy to maintain, such as those based on switchable bi- or multi-
stable circuit elements, even if some energy is required to change their value
(e.g., such as elements used in many memory devices). If constant power is
required, however, its cost could be significant, since the number of
connection elements is $M_{1}$ times larger than the number of transistors. If
connections were implemented using transistors operated as voltage-controlled
resistors, for example, the self-simulation advantages predicted above would
be up to $M_{1}$ times smaller (for energy, for speed there would be no
change). For very large transistor networks, this would not preclude enormous
computational advantages, since the self-simulation advantages found above can
very easily exceed $10^{10}$ for realistic numbers of transistors and
$M_{1}\approx 10^{3}$. Nonetheless, passive, or nearly-passive, connectivity
will be an important feature of scalable PNNs. In this regard, PAT and the PNN
concept should be useful in using coupling mechanisms that do not necessarily
exhibit ideal linearity, low noise, high-controllable-precision, or any other
number of features that would be required if devices were assumed to perfectly
obey an equation like Eqn. 44.

##  4 Supplementary discussion: when can PNNs provide technological benefits?

Drawing on the previous section and other analysis we have performed in the
preparation of this work, this section attempts to provide some heuristic
guidelines for the types of PNNs that could provide real benefits for
computation, and the conditions under which these might be achieved. This
section is not intended to be exhaustive or detailed, but merely to summarize
basic considerations for obtaining actual computational or other advantages
from PNNs.

Many tasks one could ask a PNN to do will involve a number of read-in
operations $N_{\text{RI}}$, and read-out operations $N_{\text{RO}}$. These
account for converting the data from the electronic digital domain to the PNN
domain, and from the PNN domain to the electronic digital domain. Of course,
there will also be a third set of operations, those associated with actually
performing the computation, which may be performed entirely or in part by the
physical system. If portions of the computation are performed with
conventional hardware, those operations can be counted as $N_{\text{D}}$.

To quantify the effective number of operations performed by the physical
system, $N_{\text{P,eff}}$, we can compare the PNN’s performance to an
equivalent algorithm executed on a conventional computer. For example, if a
PNN performs MNIST classification with a certain accuracy, we would say it
performs $N_{\text{P,eff}}=N_{\text{C}}$ operations, where $N_{\text{C}}$ is
the number of operations required by a digital NN that reaches the same
performance. If the PNN uses some digital operations to achieve the accuracy,
the physical portion’s effective operations count is
$N_{\text{P,eff}}=N_{\text{C}}-N_{\text{D}}$.

For computational acceleration, we ultimately care about having the total time
and energy of the PNN computation (including all parts above) be significantly
smaller than the equivalent time and energy of doing the calculation purely on
a digital computer. We therefore want the following two statements to be true.
First, to ensure an energy benefit from the PNN:

$e_{\text{RO}}N_{\text{RO}}+e_{\text{RI}}N_{\text{RI}}+e_{\text{D}}N_{\text{D}}+E_{\text{PNN}}\ll
e_{\text{D}}N_{\text{C}},$ (48)

where $e_{\text{RO}}$, $e_{\text{RI}}$, and $e_{\text{D}}$ are the per-
operation energy cost of read-out, read-in and other digital operations, and
$E_{\text{PNN}}$ is the total energy cost of executing physical operations.

Second, to ensure a speed benefit from the PNN:

$t_{\text{RO}}N_{\text{RO}}+t_{\text{RI}}N_{\text{RI}}+t_{\text{D}}N_{\text{D}}+T_{\text{PNN}}\ll
t_{\text{D}}N_{\text{C}},$ (49)

where $t_{\text{RO}}$, $t_{\text{RI}}$, and $t_{\text{D}}$ are the per-
operation time to perform read-out, read-in and other digital operations, and
$T_{\text{PNN}}$ is the total time to execute physical operations.

To a first approximation, we can assume that all digital operations, including
read-in and read-out, will take the same amount of energy and time, such that
$e_{\text{RO}}=e_{\text{RI}}=e_{\text{D}}$ and
$t_{\text{RO}}=t_{\text{RI}}=t_{\text{D}}$.

The basic idea to gain a benefit with a PNN for computation is, we expect, the
following: One wants to make sure that the task-PNN combination is
operationally bottlenecked with respect to the electronic digital computer it
is partnered with.

Operational bottleneck: A task-PNN combination is operationally bottlenecked
if the vast majority of the operations (the ‘complexity’ of the calculation
task) occurs within the PNN, with only a small amount happening at the input
(read-in), output (read-out), or on other digital operations. A task in which
we read-in $N$ numbers and find the maximum of those $N$ numbers cannot be
operationally bottlenecked, since the read-in consumes $N$ operations and the
computation only requires $\sim N$ operations. A task in which $N$ numbers are
read-in and we need to perform $\sim N^{2}$ operation-equivalents (i.e., on
the digital computer we’d have to perform $\sim N^{2}$ operations), as occurs
in an arbitrary matrix-vector operation, could in principle be operationally
bottlenecked.

With our assumptions above in place, such an operational bottleneck means:

$N_{\text{RO}}+N_{\text{RI}}+N_{\text{D}}\ll N_{\text{C}}.$ (50)

This is because the operations on the left hand side have no potential for
physical advantage: they add time and energy at the rate of any digital
computer operations. Some of the operations in the right hand side can be
executed with significant speed and/or energy advantages by the physics of the
PNN’s physical modules. Even if the energy and time required for executing the
physical system is zero, an unconventional physical processor can only provide
a benefit to the extent the above inequality is satisfied.

This is equivalent to saying that the vast majority of the computational
operations need to be performed by the physical system, since that is where
speed-ups and energy gains can be realized:

$N_{\text{P,eff}}\gg N_{\text{RO}}+N_{\text{RI}}+N_{\text{D}}.$ (51)

For example, with MNIST, a 196-dimensional image is sent into the PNN, and the
PNN produces a 10-D output, or possibly a 1-D output. If the PNN’s performance
is equivalent to a 196x10 single-layer perceptron (i.e., it gets about 90%
accuracy), then the number of equivalent-operations it performs is only 10
times higher than the number of operations required for read-in. It is thus
very hard for such a simple task to gain from a PNN: the PNN will always be
far from being operationally bottlenecked, and will fail to offer much
benefit, if any. Concretely, $N_{\text{RO}}=10$, $N_{\text{RI}}=196$, and
$N_{\text{D}}\sim 196$ (since we will almost always have to rescale data
before sending it to a physical co-processor). Meanwhile, the total number of
operations is at most 1960, so the maximum possible energy gain (and speed-up)
possible with a physical co-processor is $1960/(10+196+196)\approx 5$.

If on the other hand, the PNN is performing a much harder task, then it may be
possible for the number of operations performed physically to dwarf the read-
in and read-out operations, and any digital overhead operations. For CIFAR10,
high accuracy classification of the 3-channel, 32 by 32 images would be
equivalent to 1-100 billion operations [45]. The maximum possible PNN
advantage is then approximately $\sim 10^{9}/(32\times 32\times 3)\approx
10^{5}$.
There are several ways to overcome these challenges, which are mostly not
exclusive.
No- or few-parameter PNNs: Suppose the PNN is used as a powerful,
parameterized feature map that can be used many times within a given
calculation, or provides a complex transformation that we do not need to
control to the same level of complexity (such as if the PNN does some very
complicated nonlinear transformation, which we control by a number of
parameters much smaller than the number of equivalent digital operations that
would be required to simulate the PNN). If this feature map is useful, it can
replace a lot of digital computations that would ordinarily be required for a
task, and it may be justifiable to perform digital computations (such as input
and output linear layers) that allow the few-parameter PNN to be used in
numerous different ways within a hybrid digital-PNN computation.

This type of PNN resembles reservoir computing or kernel methods, and indeed,
the success of these methods provides some motivation as to why it would be
beneficial. However, it is important to appreciate why this is not equivalent
to reservoir computing. In reservoir computing, the physical feature map is
essentially untrained, and cannot be backpropagated through. Thus, any digital
pre-processing prior to the physical feature map cannot be efficiently
trained. In contrast, by using PAT, one may backpropagate through the physical
feature map, and thus optimize pre-processing (such as the linear input layer
used in Fig. 4K-O in the main article) to maximally utilize the physical
feature map. Using PAT, one may use sequences of physical feature maps,
possibly with digital operations between them, and train all parameters with
backpropagation.

At first glance, the optical Fourier transform performed by a lens seems to be
a good candidate for this. However, on closer inspection, it does not meet the
requirement, because one needs to read-in and read-out $N$ samples, but the
optics only provides an operations-equivalent of $N\log N$. Thus, the physical
benefit is at most around $(N\log N)/3N=(\log N)/3$, which is only significant
for extremely large $N$.

A better example is a PNN that implements a large fraction of the initial
layers of a large image-classification or speech-classification network,
generating mainly generic features [46]. This could have few trainable
parameters (even zero, if an input layer was used), but still provide a
significant benefit as long as the task would utilize the features produced by
the PNN. Here, it would still be beneficial that the PNN be differentiable, in
order to learn the input layer, or in order to train a small number of
parameters.

Parameter-in-place PNNs: If the PNN’s parameters must be read-in for each
inference, then the number of read-in operations can quickly approach the
scale of the computational operations (unless the PNN somehow only requires a
few parameters). Since read-in operations are limited by the energy and speed
of conventional processors, it is optimal that most of a PNN’s physical
parameters are stored in-place. That is, the parameters that control the PNN’s
input-output transformation may be stored in the PNN itself, rather than read-
in through the same channel as the inputs. If computations the PNN performs
are part of a typical NN inference, then the parameters do not need to be
updated as often as the inputs, and the PNN can avoid a significant read-in
cost. As was noted in the last section, it is helpful (but not essential) that
in-place parameters be fully passive, in order to minimize their energy cost.

PNNs based on imperfect physical operators: In order to minimize digital
overhead operations, one strategy would be to design physical transformations
that approximately implement specific desired operations, such as matrix-
vector multiplications, or whose controllable dynamics loosely resemble neural
networks, such as the oscillator and transistor networks discussed in the last
section. Such operations could be very approximate, with a myriad of complex
deviations from the nominal design operation. However, as physical
transformations in a PNN whose parameters (including digital pre-processing)
are going to be trained anyway, these imperfect physical operations would
likely serve their function well enough.

Physical-feedforward PNNs: In some PNNs, the outputs can be fed directly into
another PNN layer without a digital intermediary. PNNs in which the output is
directly measured (such as electronic PNNs where the output is a voltage, or
possibly in optical cavities or arrangements, where the full optical field is
fed-forward to subsequent physical processing) are in this category, as are
PNNs which have a memory of previous inputs and states (as in optical cavities
or other resonators). In this case, one can cascade PNNs over many layers,
building up computational complexity, without incurring read-in or read-out
costs except at the very beginning and end.

There are two other important scenarios where PNNs can have benefits, but that
are broader than just computational energy or speed.

Physical-data-domain-input PNNs: In some PNNs, the input data will already be
in the physical domain of the PNN, such as light entering an optical PNN. In
this case, the PNN not only avoids the read-in cost, but can additionally
perform operations that a traditional device, which would first read-in the
data to the electronic domain and then process it, could not. As in compressed
sensing single-pixel cameras using spatial light modulators [47], this could
allow for higher-speed, higher-resolution, or lower-energy costs for
performing measurements (either by improving the nature of the measurement
itself, or by performing parts of the post-processing physically). In other
words, one could use PAT to train physical ‘smart sensors’.

Physical-output PNNs: PNNs may also provide a means of producing a desired
physical quantity or functionality, in the domain of the physical system,
rather than a computation result that is ultimately utilized in a downstream
digital computer. In PNNs, the final output is normally obtained from
measuring the physical system. However, we can instead treat the final state
of the physical system as the output, and the measurement of it as merely an
assessment of that physical state. For example, a nonlinear optical PNN may be
trained to produce a desired output pulse, which would be assessed based on
measurements of that pulse (such as the power spectrum used in our optical
SHG-PNN). Alternatively, the physical functionality of the produced physical
state might be assessed, such as by measuring the effects of the produced
pulse on a chemical, the optical functionality of a complex nanostructure
[48], or the behavior of a small robot [49]. In all cases, measurements that
assess the physical state need to be performed and transferred to the digital
domain, in order to facilitate PAT, but the true output of the PNN is the
physical object and its state. Thus, in PNNs of this kind, as with the ‘smart
sensors’ described above, PAT allows the backpropagation algorithm for
gradient descent to be applied to produce a desired physical state or
functionality, in the physical domain of the PNN.

## References

  * [1] Hubara, I., Courbariaux, M., Soudry, D., El-Yaniv, R. & Bengio, Y.  Quantized neural networks: Training neural networks with low precision weights and activations.  _The Journal of Machine Learning Research_ 18, 6869–6898 (2017).
  * [2] Lillicrap, T. P., Cownden, D., Tweed, D. B. & Akerman, C. J.  Random synaptic feedback weights support error backpropagation for deep learning.  _Nature Communications_ 7, 1–10 (2016).
  * [3] Nielsen, M.  _Neural Networks and Deep Learning_ (Determination Press, 2015).
  * [4] Kingma, D. P. & Ba, J.  Adam: A method for stochastic optimization.  In Bengio, Y. & LeCun, Y. (eds.) _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ (2015).
  * [5] Zeiler, M. D.  Adadelta: an adaptive learning rate method.  _arXiv:1212.5701_ (2012).
  * [6] Liu, W. et al.  Programmable controlled mode-locked fiber laser using a digital micromirror device.  _Optics Letters_ 42, 1923–1926 (2017).
  * [7] Matthès, M. W., del Hougne, P., de Rosny, J., Lerosey, G. & Popoff, S. M.  Optical complex media as universal reconfigurable linear operators.  _Optica_ 6, 465–472 (2019).
  * [8] Popoff, S. M. & Matthès, M. W.  ALP4lib: A Python wrapper for the Vialux ALP-4 controller suite to control DMDs, https://doi.org/10.5281/zenodo.4076193 (2020).
  * [9] Akiba, T., Sano, S., Yanase, T., Ohta, T. & Koyama, M.  Optuna: A next-generation hyperparameter optimization framework.  In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, 2623–2631 (2019).
  * [10] Torrejon, J. et al.  Neuromorphic computing with nanoscale spintronic oscillators.  _Nature_ 547, 428–431 (2017).
  * [11] Veit, A., Wilber, M. & Belongie, S.  Residual networks behave like ensembles of relatively shallow networks.  _arXiv:1605.06431_ (2016).
  * [12] Lin, H. W., Tegmark, M. & Rolnick, D.  Why does deep and cheap learning work so well?  _Journal of Statistical Physics_ 168, 1223–1247 (2017).
  * [13] History of supercomputing, https://en.wikipedia.org/wiki/History_of_supercomputing (2020).
  * [14] Green500 November 2020, https://www.top500.org/lists/green500/2020/11/ (2020).
  * [15] Geforce 20 series - Wikipedia, https://en.wikipedia.org/wiki/GeForce_20_series (2021).
  * [16] Agrawal, G. P.  _Nonlinear Fiber Optics_ (Elsevier, 2000).
  * [17] Carlson, D. R. et al.  Ultrafast electro-optic light with subcycle control.  _Science_ 361, 1358–1363 (2018).
  * [18] Jankowski, M. et al.  Ultrabroadband nonlinear optics in nanophotonic periodically poled lithium niobate waveguides.  _Optica_ 7, 40–46 (2020).
  * [19] Xu, B., Zhou, Y. & Chiu, Y.  A 23mW 24GS/s 6b time-interleaved hybrid two-step ADC in 28nm CMOS.  In _2016 IEEE Symposium on VLSI Circuits (VLSI-Circuits)_ , 1–2 (IEEE, 2016).
  * [20] Sun, C. et al.  Single-chip microprocessor that communicates directly using light.  _Nature_ 528, 534–538 (2015).
  * [21] Atabaki, A. H. et al.  Integrating photonics with silicon nanoelectronics for the next generation of systems on a chip.  _Nature_ 556, 349–354 (2018).
  * [22] Wang, C. et al.  Integrated lithium niobate electro-optic modulators operating at CMOS-compatible voltages.  _Nature_ 562, 101–104 (2018).
  * [23] Koeber, S. et al.  Femtojoule electro-optic modulation using a silicon–organic hybrid device.  _Light: Science & Applications_ 4, e255–e255 (2015).
  * [24] Srinivasan, S. A. et al.  56 Gb/s germanium waveguide electro-absorption modulator.  _Journal of Lightwave Technology_ 34, 419–424 (2015).
  * [25] Nozaki, K. et al.  Photonic-crystal nano-photodetector with ultrasmall capacitance for on-chip light-to-voltage conversion without an amplifier.  _Optica_ 3, 483–492 (2016).
  * [26] Tang, L. et al.  Nanometre-scale germanium photodetector enhanced by a near-infrared dipole antenna.  _Nature Photonics_ 2, 226–229 (2008).
  * [27] Rafayelyan, M., Dong, J., Tan, Y., Krzakala, F. & Gigan, S.  Large-scale optical reservoir computing for spatiotemporal chaotic systems prediction.  _Physical Review X_ 10, 041037 (2020).
  * [28] Teğin, U., Yıldırım, M., Oğuz, İ., Moser, C. & Psaltis, D.  Scalable optical learning operator.  _arXiv:2012.12404_ (2020).
  * [29] Popoff, S. M. et al.  Measuring the transmission matrix in optics: an approach to the study and control of light propagation in disordered media.  _Physical Review Letters_ 104, 100601 (2010).
  * [30] Xiong, W. et al.  Principal modes in multimode fibers: exploring the crossover from weak to strong mode coupling.  _Optics Express_ 25, 2709–2724 (2017).
  * [31] Ohana, R. et al.  Kernel computations from large-scale random features obtained by optical processing units.  In _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 9294–9298 (IEEE, 2020).
  * [32] Kasim, M. et al.  Up to two billion times acceleration of scientific simulations with deep neural architecture search.  _arXiv:2001.08055_ (2020).
  * [33] Poletti, F. & Horak, P.  Description of ultrashort pulse propagation in multimode optical fibers.  _JOSA B_ 25, 1645–1654 (2008).
  * [34] Lægsgaard, J.  Efficient simulation of multimodal nonlinear propagation in step-index fibers.  _JOSA B_ 34, 2266–2273 (2017).
  * [35] Wright, L. G. et al.  Multimode nonlinear fiber optics: massively parallel numerical solver, tutorial, and outlook.  _IEEE Journal of Selected Topics in Quantum Electronics_ 24, 1–16 (2017).
  * [36] Sato, H. et al.  Properties of magnetic tunnel junctions with a MgO/CoFeB/Ta/CoFeB/MgO recording structure down to junction diameter of 11 nm.  _Applied Physics Letters_ 105, 062403 (2014).
  * [37] He, K., Zhang, X., Ren, S. & Sun, J.  Deep residual learning for image recognition.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 770–778 (2016).
  * [38] Chen, R. T., Rubanova, Y., Bettencourt, J. & Duvenaud, D.  Neural ordinary differential equations.  _arXiv:1806.07366_ (2018).
  * [39] Günther, M., Feldmann, U. & ter Maten, J.  Modelling and discretization of circuit problems.  _Handbook of Numerical Analysis_ 13, 523–659 (2005).
  * [40] Rewieński, M.  A perspective on fast-spice simulation technology.  In _Simulation and Verification of Electronic and Biological Systems_ , 23–42 (Springer, 2011).
  * [41] Benk, J., Denk, G. & Waldherr, K.  A holistic fast and parallel approach for accurate transient simulations of analog circuits.  _Journal of Mathematics in Industry_ 7, 1–19 (2017).
  * [42] Sun, Y., Agostini, N. B., Dong, S. & Kaeli, D.  Summarizing CPU and GPU design trends with product data.  _arXiv:1911.11313_ (2019).
  * [43] Enz, C. C., Krummenacher, F. & Vittoz, E. A.  An analytical MOS transistor model valid in all regions of operation and dedicated to low-voltage and low-current applications.  _Analog integrated circuits and signal processing_ 8, 83–114 (1995).
  * [44] Enz, C. C.  An MOS transistor model for RF IC design valid in all regions of operation.  _IEEE Transactions on Microwave Theory and Techniques_ 50, 342–359 (2002).
  * [45] Hasanpour, S. H., Rouhani, M., Fayyaz, M. & Sabokrou, M.  Lets keep it simple, using simple architectures to outperform deeper and more complex architectures.  _arXiv:1608.06037_ (2016).
  * [46] Sharif Razavian, A., Azizpour, H., Sullivan, J. & Carlsson, S.  CNN features off-the-shelf: an astounding baseline for recognition.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops_ , 806–813 (2014).
  * [47] Duarte, M. F. et al.  Single-pixel imaging via compressive sampling.  _IEEE signal processing magazine_ 25, 83–91 (2008).
  * [48] Peurifoy, J. et al.  Nanophotonic particle simulation and inverse design using artificial neural networks.  _Science Advances_ 4, eaar4206 (2018).
  * [49] Howison, T., Hauser, S., Hughes, J. & Iida, F.  Reality-assisted evolution of soft robots through large-scale physical experimentation: a review.  _Artificial Life_ 26, 484–506 (2021).