Quantitative Approximation for Neural Operators in Nonlinear Parabolic Equations

## Abstract

Abstract Neural operators serve as universal approximators for general continuous operators. In this paper, we derive the approximation rate of solution operators for the nonlinear parabolic partial differential equations (PDEs), contributing to the quantitative approximation theorem for solution operators of nonlinear PDEs. Our results show that neural operators can efficiently approximate these solution operators without the exponential growth in model complexity, thus strengthening the theoretical foundation of neural operators. A key insight in our proof is to transfer PDEs into the corresponding integral equations via Duahamel’s principle, and to leverage the similarity between neural operators and Picard’s iteration—a classical algorithm for solving PDEs. This approach is potentially generalizable beyond parabolic PDEs to a range of other equations, including the Navier-Stokes equation, nonlinear Schrödinger equations and nonlinear wave equations, which can be solved by Picard’s iteration.

## 1 Introduction

Neural operators have gained significant attention in deep learning as an extension of traditional neural networks. While conventional neural networks are designed to learn mappings between finite-dimensional spaces, neural operators extend this capability by learning mappings between infinite-dimensional function spaces.
A key application of neural operators is in constructing surrogate models of solvers for partial differential equations (PDEs) by learning their solution operators. Traditional PDE solvers often require substantial computational resources and time, especially when addressing problems with high dimensionality, nonlinearity, or complex boundary shape. In contrast, once trained, neural operators serve as surrogate models, providing significantly faster inference compared to traditional numerical solvers.
Neural operators are recognized as universal approximators for general continuous operators. However, the theoretical understanding of their approximation capabilities, particularly as solvers for PDEs, is not yet fully developed.

This paper focuses on whether neural operators are suitable for approximating solution operators of PDEs, and which neural operator architectures might be effective for this purpose. The idea of this paper is to define neural operators by aligning them with Picard’s iteration, a classical method for solving PDEs. Specifically, by associating each forward pass of the neural operator’s layers with a Picard’s iteration step, we hypothesize that increasing the number of layers would naturally lead to an approximate solution of the PDE. We constructively prove a quantitative approximation theorem for solution operators of PDEs based on this idea. This theorem shows that, by appropriately selecting the basis functions within the neural operator, it is possible to avoid the exponential growth in model complexity (also called “the curse of parametric complexity”) that often arises in general operator approximation, thereby providing a theoretical justification for the effectiveness of neural operators as PDE solvers.

### 1.1 Related works

Neural operators were introduced by [] as one of the operator learning methods, such as DeepONet [, ] and PCA-Net []. Various architectures have been proposed, including Graph Neural Operators [], Fourier Neural Operators [], Wavelet Neural Operators [, ], Spherical Fourier Neural Operators [], and Laplace Neural Operators []. These architectures have demonstrated empirical success as the surrogate models of simulators across a wide range of PDEs, as benchmarked in [].
In the case of parabolic PDEs, which are the focus of this paper, neural operators have also shown promising results, for example,
the Burgers, Darcy, Navier-Stokes equations [],
the KPP-Fisher equation [], the Allen-Cahn equation [, ], and the Nagumo equation [].

Universal approximation theorems of operator learning for general operators were established for neural operators [, , , ], DeepOnet [, ], and PCA-net []. This indicates that these learning methods possess the capabilities to approximate a wide range of operators. However, operator learning for general operators suffers from “the curse of parametric complexity”, where the number of learnable parameters exponentially grows as the desired approximation accuracy increases [].

A common approach to mitigating the curse of parametric complexity is to restrict general operators to the solution operators of PDEs. Recently, several quantitative approximation theorems have been established for the solution operators of specific PDEs without experiencing exponential growth in model complexity. For instance, [, ] developed quantitative approximation theorems for the Darcy and Navier-Stokes equations using Fourier neural operators and PCA-Net, respectively. Additionally, [] addressed the Hamilton-Jacobi equations with Hamilton-Jacobi neural operators. Further studies, such as [, , ], investigated quantitative approximation theorems using DeepONet for a range of PDEs, including elliptic, parabolic, and hyperbolic equations, while [] focused on advection-diffusion equations.

In line with these researches, the present paper also concentrates on restricting the learning operator to the solution operators of specific PDEs, namely parabolic PDEs. However, unlike previous studies, this work leverages the similarity with the Picard’s iteration in the framework of relatively general neural operator. This approach is potentially generalizable beyond parabolic PDEs to a range of other equations, including the Navier-Stokes equation, nonlinear dispersive equations, and nonlinear hyperbolic equations, which are solvable by Picard’s iteration.
See [] for other directions on the quantitative approximation of operator learning.

### 1.2 Our results and contributions

In this paper, we present a quantitative approximation theorem for the solution operator of nonlinear parabolic PDEs using neural operators. Notably, in Theorem 1, we show that for any given accuracy, the depth and the number of neurons of the neural operators do not grow exponentially.

The proof relies on Banach’s fixed point theorem. Under appropriate conditions, the solution to nonlinear PDEs can be expressed as the fixed point of a contraction mapping, which can be implemented through Picard’s iteration. This contraction mapping corresponds to an integral operator whose kernel is the Green function associated with the linear equation. By expanding the Green function in a certain basis and truncating the expansion, we approximate the contraction mapping using a layer of neural operator. In this framework, the forward propagation through the layers is interpreted as steps in Picard’s iteration.
Our approach does not heavily rely on the universality of neural networks. More specifically, we utilize their universality only to approximate the nonlinearity, which is represented as a one-dimensional nonlinear function, and hence, the approximation rates obtained in Theorem 1 correspond to the one-dimensional approximation rates of universality.
As a result, our approach demonstrates that exponential growth in model complexity of our neural operators can be avoided.

### 1.3 Notation

We introduce the notation often used in this paper.

- •
For $r\in[1,\infty]$, we write the Hölder conjugate of $r$ as $r^{\prime}$: $1/r+1/r^{\prime}=1$ if $r\in(1,\infty)$;
$r^{\prime}=\infty$ if $r=1$;
$r^{\prime}=1$ if $r=\infty$.
- •
Let $X$ be a set.
For an operator $A:X\to X$ and $k\in\mathbb{N}$, we denote by $A^{[k]}$ the $k$ times compositions of $A$ (or the $k$ times products of $A$): $A^{[0]}$ means the identity operator on $X$ and
$A^{[k]}:=\underbrace{A\circ\cdots\circ A}_{\text{$k$ times}}.$
- •
Let $X$ and $Y$ be normed spaces with norm $\|\cdot\|_{X}$ and $\|\cdot\|_{Y}$, respectively.
For a linear operator $A:X\to Y$,
we denote by $\|A\|_{X\to Y}$ the operator norm of $A$:
$\|A\|_{X\to Y}:=\sup_{\|f\|_{X}=1}\|Af\|_{Y}.$
- •
Let $X$ be a normed space with norm $\|\cdot\|_{X}$. We denote by $B_{X}(R)$ the closed ball in $X$ with center 0 0 and radius $R>0$:
$B_{X}(R):=\{f\in X:\|f\|_{X}\leq R\}.$
- •
For $q\in[1,\infty]$,
the Lebesgue space $L^{q}(D)$ is defined by the set of all measurable functions $f=f(x)$ on $D$ such that
$\|f\|_{L^{q}}:=\begin{dcases}\left(\int_{D}|f(x)|^{q}\,dx\right)^{\frac{1}{q}}
<\infty\quad&\text{if }q\not=\infty,\\
\underset{x\in D}{\operatorname{ess\ sup}}|f(x)|<\infty\quad&\text{if }q=
\infty.\end{dcases}$
For $r,s\in[1,\infty]$, the space $L^{r}(0,T;L^{s}(D))$ is defined by the set of all measurable functions $f=f(t,x)$ on $(0,T)\times D$ such that
$\|f\|_{L^{r}(0,T;L^{s})}:=\left\{\int_{0}^{T}\left(\int_{D}|f(t,x)|^{s}\,dx
\right)^{\frac{r}{s}}\,dt\right\}^{\frac{1}{r}}<\infty$
(with the usual modifications for $r=\infty$ or $s=\infty$).
- •
For $r,s\in[1,\infty]$,
the notation $\langle\cdot,\cdot\rangle$ means the dual pair of $L^{s^{\prime}}(D)$ and $L^{s}(D)$ or $L^{r^{\prime}}(0,T;L^{s^{\prime}}(D))$ and $L^{r}(0,T;L^{s}(D))$:
$\langle u,v\rangle:=\int_{D}u(x)v(x)\,dx\quad\text{or}\quad\langle u,v\rangle:
=\int_{0}^{T}\int_{D}u(t,x)v(t,x)\,dxdt,$
respectively.

## 2 Local well-posedness for nonlinear parabolic equations

To begin with, we describe the problem setting of PDEs addressed in this paper.
Let $D$ be a bounded domain in $\mathbb{R}^{d}$ with $d\in\mathbb{N}$.
We consider
the Cauchy problem for the following nonlinear parabolic PDEs:

$$ $\begin{cases}\partial_{t}u+\mathcal{L}u=F(u)\quad&\text{in }(0,T)\times D,\\ u(0)=u_{0}\quad&\text{in }D,\end{cases}$ (P) $$

where $T>0$, $\partial_{t}:=\partial/\partial t$ denotes the time derivative,
$\mathcal{L}$ is a certain operator (e.g. $\mathcal{L}=-\Delta:=-\sum_{j=1}^{d}\partial^{2}/\partial x_{j}^{2}$ is the Laplacian),
$F:\mathbb{R}\to\mathbb{R}$ is a nonlinearity,
$u:(0,T)\times D\to\mathbb{R}$ is a solution to (P) (unknown function), and $u_{0}:D\to\mathbb{R}$ is an initial data (prescribed function).
It is important to note that boundary conditions are contained in the operator $\mathcal{L}$.
In this sense, the problem (P) can be viewed as an abstract initial boundary value problem on $D$. See Appendix A for examples of $\mathcal{L}$.

As it is a standard practice, we study
the problem (P) via the integral formulation

$$ $u(t)=S_{\mathcal{L}}(t)u_{0}+\int_{0}^{t}S_{\mathcal{L}}(t-\tau)F(u(\tau))\,d\tau,$ (P’) $$

where $\{S_{\mathcal{L}}(t)\}_{t\geq 0}$ is the semigroup
generated by $\mathcal{L}$ (i.e. the solution operators of the linear equation $\partial_{t}u+\mathcal{L}u=0$).
The second term in the right hand side of (P’) is
commonly referred to as
the Duhamel’s integral.
Under suitable conditions on $\mathcal{L}$ and $F$, the problems (P) and (P’) are equivalent if the function $u$ is sufficiently smooth.

The aim of this section is to state the results on local well-posedness (LWP) for (P’). Here, the LWP means the
existence of local in time solution, the uniqueness of the solution, and the continuous dependence on initial data.
Its proof is based on the fixed point argument (or also called the contraction mapping argument).
These results are fundamental to study our neural operator in Sections 3 and 4.
In particular, Propositions 1 and 2 below serve as guidelines for setting function spaces as the domain and range of neural operators in Definition 1 (Section 3) and for determining the norm to measure the error in Theorem 1 (Section 4).

In this paper we impose the following assumptions on $\mathcal{L}$ and $F$.

###### Assumption 1 .

For any $1\leq q_{1}\leq q_{2}\leq\infty$,
there exists a constant $C_{\mathcal{L}}>0$ such that

$$ $\|S_{\mathcal{L}}(t)\|_{L^{q_{1}}\to L^{q_{2}}}\leq C_{\mathcal{L}}t^{-\nu( \frac{1}{q_{1}}-\frac{1}{q_{2}})},\quad t\in(0,1],$ (1) $$

for some $\nu>0$.

###### Assumption 2 .

$F\in C^{1}(\mathbb{R};\mathbb{R})$ satisfies $F(0)=0$ and

$$ $|F(z_{1})-F(z_{2})|\leq C_{F}\max_{i=1,2}|z_{i}|^{p-1}|z_{1}-z_{2}|,$ (2) $$

for any $z_{1},z_{2}\in\mathbb{R}$ and for some $p>1$ and $C_{F}>0$.

###### Remark 1 .

The range $(0,1]$ of $t$ in (1) can be generalized to $(0,T_{\mathcal{L}}]$, but it is assumed here to be $(0,1]$ for simplicity. This generalization is not essential, as the existence time $T$ of the solution is sufficiently small in the fixed point argument later. Long time solutions are achieved by repeatedly using the solution operator of (P’) constructed in the fixed point argument (see also Subsection 4.2).

###### Remark 2 .

Typical examples of $\mathcal{L}$ are the Laplacian with the Dirichlet, Neumann, or Robin boundary condition, the Schrödinger operator, the elliptic operato,r and the higher-order Laplacian. See Appendix A for the details.
On the other hand, typical examples of the nonlinearity $F$ are $F(u)=\pm|u|^{p-1}u,\pm|u|^{p}$ (which can be regarded as the main term of the Taylor expansion of a more general nonlinearity $F$ if $F$ is smooth in some extent).

Under the above assumptions,
the problem (P’) is local well-posed, where, as a solution space, we use the space $L^{r}(0,T;L^{s}(D))$ with
the parameters $r,s$ satisfying

$$ $r,s\in[p,\infty]\quad\text{and}\quad\frac{\nu}{s}+\frac{1}{r}<\frac{1}{p-1}.$ (3) $$

More precisely, we have the following result on LWP.

###### Proposition 1 .

Assume that $r,s$ satisfy
(3). Then,
for any $u_{0}\in L^{\infty}(D)$, there exist a time $T=T(u_{0})\in(0,1]$ and a unique solution $u\in L^{r}(0,T;L^{s}(D))$ to (P’). Moreover, for any $u_{0},v_{0}\in L^{\infty}(D)$, the solutions $u$ and $v$ to (P’) with $u(0)=u_{0}$ and $v(0)=v_{0}$ satisfy the continuous dependence on initial data:
There exists a constant $C>0$ such that

$$ $\|u-v\|_{L^{r}(0,T^{\prime};L^{s})}\leq C\|u_{0}-v_{0}\|_{L^{\infty}},$ $$

where $T^{\prime}<\min\{T(u_{0}),T(v_{0})\}$.

The proof is based on the fixed point argument. Given an initial data $u_{0}\in L^{\infty}(D)$ and $T,M>0$, we define the map $\Phi=\Phi_{u_{0}}$ by

$$ $\Phi[u](t):=S_{\mathcal{L}}(t)u_{0}+\int_{0}^{t}S_{\mathcal{L}}(t-\tau)F(u( \tau))\,d\tau$ (4) $$

for $t\in[0,T]$ and the complete metric space $X:=B_{L^{r}(0,T;L^{s}(D))}(M)$
equipped with the metric

$$ $\mathrm{d}(u,v):=\|u-v\|_{L^{r}(0,T;L^{s})}.$ $$

Let $R,M>0$ be arbitrarily fixed, and let $T>0$ (which is taken sufficiently small later).
Then, under Assumptions 1 and 2, it can be proved that
for any $u_{0}\in B_{L^{\infty}}(R)$, the map $\Phi:X\to X$ is $\delta$-contractive with a contraction rate $\delta\in(0,1)$, i.e.

$$ $\mathrm{d}(\Phi[u],\Phi[v])\leq\delta\mathrm{d}(u,v)\quad\text{for any $u,v\in X $},$ $$

where $T$ is taken small enough to depend on $R$, $M$ and $\delta$ (not to depend on $u_{0}$ itself).
Therefore, Banach’s fixed point theorem allows us to prove that there exists uniquely
a function $u\in X$ such that $u=\Phi[u]$ (a fixed point).
This function $u$ is precisely the solution of (P’) with the initial data $u(0)=u_{0}$. Thus, Proposition 1 is shown. See Appendix B for more details of the proof.

This proof by the fixed point argument guarantees the following result (see e.g. []).

###### Proposition 2 .

Assume that $r,s$ satisfy
(3). Then, for any $R,M>0$ and for any $\delta\in(0,1)$, there exists a time $T\in(0,1]$, depending on $R,M$ and $\delta$,
such that
the following statements hold:

- (i)
There exists a unique solution operator $\Gamma^{+}:B_{L^{\infty}}(R)\to B_{L^{r}(0,T;L^{s})}(M)$
such that
$\Gamma^{+}(u_{0})=u$
for any $u_{0}\in B_{L^{\infty}}(R)$, where $u$ is the solution to (P’) with $u(0)=u_{0}$ given in Proposition 1.
- (ii)
Given $u_{0}\in B_{L^{\infty}}(R)$, define Picard’s iteration by
$\begin{dcases}u^{(1)}:=S_{\mathcal{L}}(t)u_{0},\\
u^{(\ell+1)}:=\Phi[u^{(\ell)}]=u^{(1)}+\int_{0}^{t}S_{\mathcal{L}}(t-\tau)F(u^
{(\ell)}(\tau))\,d\tau,\quad\ell=1,2,\cdots,\end{dcases}$
(5)
that is,
$u^{(\ell)}:=\Phi^{[\ell]}[0]=\underbrace{\Phi\circ\cdots\circ\Phi}_{\text{$
\ell$ times}}[0],\quad\ell=1,2,\cdots,$
where $\Phi:X\to X$ is a $\delta$-contraction mapping defined by (4).
Then $u^{(\ell)}\to u$ in $L^{r}(0,T;L^{s}(D))$ as $\ell\to\infty$ and
$\mathrm{d}(u^{(\ell)},u)\leq\frac{\delta^{\ell}}{1-\delta}\mathrm{d}(u^{(1)},0
),\quad\ell=1,2,\cdots.$

###### Assumption 1 .

###### Assumption 2 .

###### Remark 1 .

###### Remark 2 .

###### Proposition 1 .

###### Proposition 2 .

## 3 Neural operator for nonlinear parabolic equations

In this section, we aim to construct neural operators $\Gamma$
that serve as accurate approximation models of the solution operators $\Gamma^{+}$
for nonlinear parabolic PDEs (P).
We start by explaining our idea in rough form.
Our idea is inspired by the fixed point argument and Picard’s iteration.
By Section 2, for any $u_{0}\in B_{L^{\infty}}(R)$, the solution $u$ to (P’) on $[0,T]\times D$ can be obtained through the Picard’s iteration (5) under appropriate settings.
If the semigroup $S_{\mathcal{L}}(t)$ has an integral kernel $G=G(t,x,y)$, which is the Green function $G$ of the linear equation $\partial_{t}u+\mathcal{L}u=0$, then we can write

$$ $S_{\mathcal{L}}(t)u_{0}(x)=\int_{D}G(t,x,y)u_{0}(y)\,dy,$ $$

$$ $\int_{0}^{t}S_{\mathcal{L}}(t-\tau)F(u^{(\ell)}(\tau,x))\,d\tau=\int_{0}^{t} \int_{D}G(t-\tau,x,y)F(u^{(\ell)}(\tau,y))\,dyd\tau.$ $$

In addition, we suppose that the kernel $G$ has an expansion

$$ $\begin{split}G(t-\tau,x,y)&=\sum_{m,n\in\Lambda}c_{m,n}\psi_{m}(\tau,y)\varphi _{n}(t,x)=\lim_{N\to\infty}\sum_{m,n\in\Lambda_{N}}c_{m,n}\psi_{m}(\tau,y) \varphi_{n}(t,x),\end{split}$ $$

for $0\leq\tau<t\leq T$ and $x,y\in D$,
where $\Lambda$ is an index set that is either finite or countably infinite, and $\Lambda_{N}$ is a subset of $\Lambda$ with its cardinality $|\Lambda_{N}|=N\in\mathbb{N}$ and the monotonicity $\Lambda_{N}\subset\Lambda_{N^{\prime}}$ for any $N\leq N^{\prime}$.
We write the partial sum as

$$ $G_{N}(t-\tau,x,y):=\sum_{m,n\in\Lambda_{N}}c_{m,n}\psi_{m}(\tau,y)\varphi_{n}( t,x)$ $$

for $0\leq\tau<t\leq T$ and $x,y\in D$.
We define $\Phi_{N}$ by

$$ $\displaystyle\Phi_{N}[u](t,x):=$ $\displaystyle\int_{D}G_{N}(t,x,y)u_{0}(y)\,dy+\int_{0}^{t}\int_{D}G_{N}(t-\tau ,x,y)F(u(\tau,y))\,dyd\tau$ $\displaystyle=$ $\displaystyle\sum_{m,n\in\Lambda_{N}}c_{m,n}\langle\psi_{m}(0,\cdot),u_{0} \rangle\varphi_{n}(t,x)+\sum_{m,n\in\Lambda_{N}}c_{m,n}\langle\psi_{m},F(u) \rangle\varphi_{n}(t,x)$ $\displaystyle=$ $\displaystyle\sum_{m,n\in\Lambda_{N}}c_{m,n}\left(\langle\psi_{m}(0,\cdot),u_{ 0}\rangle+\langle\psi_{m},F(u)\rangle\right)\varphi_{n}(t,x),$ $$

and moreover, we define an approximate Picard’s iteration by

$$ $\begin{dcases}\hat{u}^{(1)}:=\int_{D}G_{N}(t,x,y)u_{0}(y)\,dy=\sum_{m,n\in \Lambda_{N}}c_{m,n}\langle\psi_{m}(0,\cdot),u_{0}\rangle\varphi_{n}(t,x),\\ \hat{u}^{(\ell+1)}:=\Phi_{N}[\hat{u}^{(\ell)}],\quad\ell=1,2,\cdots.\end{dcases}$ $$

Then, for sufficiently large $N$, we expect:

- 1.
$\Phi_{N}\approx\Phi$ and $\Phi_{N}$ is also contractive on $X_{N}:=B_{L^{r}(0,T_{N};L^{s})}(M)$ for some $T_{N}>0$.
- 2.
There exists a fixed point $\hat{u}\in X_{N}$ of $\Phi_{N}$ such that
$\hat{u}^{(\ell)}\to\hat{u}$ as $\ell\to\infty$.
- 3.
$\hat{u}\approx u$ on $[0,\min\{T,T_{N}\}]\times D$ for any $u_{0}\in B_{L^{\infty}}(R)$.
- 4.
Define $\Gamma:B_{L^{\infty}}(R)\to X_{N}$ by $\Gamma(u_{0}):=\hat{u}^{(L)}$ for $L\in\mathbb{N}$. Then
$\Gamma\approx\Gamma^{+}$ as $N,L\to\infty$.

The above $\Gamma$ is the prototype of our neural operator, where $N$ corresponds to the rank and $L$ to the layer depth.
In other words, in our neural operator, the Picard’s iteration step corresponds to the forward propagation in the layer direction of the neural operator, which converges to the fixed point (i.e. the solution to (P’)) as the layers get deeper for sufficiently large rank $N$.

The above is only a rough idea.
The precise definition of our neural operator is the following:

###### Definition 1 (Neural operator) .

Let $T>0$ and $r,s\in[1,\infty]$.
Let $\varphi:=\{\varphi_{n}\}_{n}$ and $\psi:=\{\psi_{m}\}_{m}$ be families of functions in
$L^{r}(0,T;L^{s}(D))$ and $L^{r^{\prime}}(0,T;L^{s^{\prime}}(D))$, respectively.
We define a neural operator $\Gamma:L^{\infty}(D)\to L^{r}(0,T;L^{s}(D))$ by

$$ $\Gamma:L^{\infty}(D)\to L^{r}(0,T;L^{s}(D)):u_{0}\mapsto\hat{u}^{(L+1)}.$ $$

Here, the output function $\hat{u}^{(L+1)}$ is given by the following steps:

1. (Input layer)
$\hat{u}^{(1)}=(\hat{u}^{(1)}_{1},\hat{u}^{(1)}_{2},\ldots,\hat{u}^{(1)}_{d_{1}})$ is given by

$$ $\hat{u}^{(1)}(t,x):=(K^{(0)}_{N}u_{0})(t,x)+b^{(0)}_{N}(t,x).$ $$

Here, $K^{(0)}_{N}:L^{\infty}(D)\to L^{r}(0,T;L^{s}(D))^{d_{1}}$ and $b^{(0)}_{N}\in L^{r}(0,T;L^{s}(D))^{d_{1}}$ are defined by

$$ $\displaystyle(K^{(0)}_{N}u_{0})(t,x):=\sum_{m,n\in\Lambda_{N}}C_{n,m}^{(0)} \langle\psi_{m}(0,\cdot),u_{0}\rangle\varphi_{n}(t,x)\quad\text{with $C^{(0)}_ {n,m}\in\mathbb{R}^{d_{1}\times 1}$},$ $\displaystyle b^{(0)}_{N}(t,x):=\sum_{n\in\Lambda_{N}}b_{N}^{(0)}\varphi_{n}(t ,x)\quad\text{with $b^{(0)}_{N}\in\mathbb{R}^{d_{1}}$}.$ $$

2. (Hidden layers)
For $2\leq\ell\leq L$, $\hat{u}^{(\ell)}=(\hat{u}^{(\ell)}_{1},\hat{u}^{(\ell)}_{2},\ldots,\hat{u}^{(
\ell)}_{d_{\ell}})$ are iteratively given by

$$ $\hat{u}^{(\ell+1)}(t,x)=\sigma\left(W^{(\ell)}\hat{u}^{(\ell)}(t,x)+(K^{(\ell) }_{N}\hat{u}^{(\ell)})(t,x)+b^{(\ell)}\right),\quad 1\leq\ell\leq L-1.$ $$

3. (Output layer)
$\hat{u}^{(L+1)}$ is given by

$$ $\hat{u}^{(L+1)}(t,x)=W^{(L)}\hat{u}^{(L)}(t,x)+(K^{(L)}_{N}\hat{u}^{(L)})(t,x) +b^{(L)}.$ $$

Here, $\sigma:\mathbb{R}\to\mathbb{R}$ is a nonlinear activation operating element-wise, and $W^{(\ell)}\in\mathbb{R}^{d_{\ell+1}\times d_{\ell+1}}$ is a weight matrix of the $\ell$-th hidden layer, and $b^{(\ell)}\in\mathbb{R}^{d_{\ell}}$ is a bias vector, and $K^{(\ell)}_{N}:L^{r}(0,T;L^{s}(D))^{d_{\ell}}\to L^{r}(0,T;L^{s}(D))^{d_{\ell+
1}}$
is defined by

$$ $\displaystyle(K^{(\ell)}_{N}u)(t,x):=\sum_{m,n\in\Lambda_{N}}C_{n,m}^{(\ell)} \langle\psi_{m},u\rangle\varphi_{n}(t,x)\quad\text{with $C^{(\ell)}_{n,m}\in \mathbb{R}^{d_{\ell+1}\times d_{\ell}}$},$ $$

where we use the notation

$$ $\langle\psi_{m},u\rangle:=\left(\langle\psi_{m},u_{1}\rangle,\ldots,\langle \psi_{m},u_{d_{\ell}}\rangle\right)\in\mathbb{R}^{d_{\ell}},$ $$

for $u=(u_{1},\dots,u_{d_{\ell}})\in L^{r}(0,T;L^{s}(D))^{d_{\ell}}$.
Note that $d_{L+1}=1$.

We denote by $\mathcal{NO}^{L,H,\sigma}_{N,\varphi,\psi}$ the class of neural operators defined as above, with the depth $L$, the number of neurons $H=\sum_{\ell=1}^{L}d_{\ell}$, the rank $N$, the activation function $\sigma$, and the families of functions $\varphi,\psi$.

Note that $\varphi$ and $\psi$ are hyperparameters.
The examples are the Fourier basis, wavelet basis, orthogonal polynomial, spherical harmonics, and eigenfunctions of $\mathcal{L}$.
If we select $\varphi$ and $\psi$ as the Fourier
basis and wavelet basis,
the network architectures correspond to the Fourier neural operator (FNO) in [] and the wavelet neural operator (WNO) in [], respectively.
See also Appendix C for FNOs and WNOs in some simple cases.

###### Definition 1 (Neural operator) .

## 4 Quantitative approximation theorem

In this section we prove a quantitative approximation theorem for our neural operator $\Gamma$ (Definition 1). For this purpose, we assume that $\mathcal{L}$ and $F$ satisfy
Assumptions 1 and 2, respectively, and the parameters $r,s$ satisfy the condition (3).
In addition,
we assume the following:

###### Assumption 3 .

Let $\Lambda$ be an index set that is either finite or countably infinite and $\Lambda_{N}$ a subset of $\Lambda$ with its cardinality $|\Lambda_{N}|=N\in\mathbb{N}$ and the monotonicity $\Lambda_{N}\subset\Lambda_{N^{\prime}}$ for any $N\leq N^{\prime}$.
For the Green function $G$ of the linear equation $\partial_{t}u+\mathcal{L}u=0$, there exist the families of functions $\varphi:=\{\varphi_{n}\}_{n\in\Lambda}$ and $\psi:=\{\psi_{m}\}_{m\in\Lambda}$ in
$L^{r}(0,T;L^{s}(D))$ and $L^{r^{\prime}}(0,T;L^{s^{\prime}}(D))$, respectively, such that $G$ has the following expansions:

$$ $G(t-\tau,x,y)=\sum_{m,n\in\Lambda}c_{n,m}\psi_{m}(\tau,y)\varphi_{n}(t,x), \quad 0<\tau<t<T,\ x,y\in D,$ (6) $$

$$ $G(t,x,y)=\sum_{m,n\in\Lambda}c_{n,m}\psi_{m}(0,y)\varphi_{n}(t,x),\quad 0<t<T, \ x,y\in D$ (7) $$

in the sense that

$$ $\displaystyle C_{G}(N):=\left\|\left\|G(t-\tau,x,y)-G_{N}(t-\tau,x,y)\right\|_ {L^{r^{\prime}}_{\tau}(0,t;L^{s^{\prime}}_{y})}\right\|_{L^{r}_{t}(0,T;L^{s}_{ x})}\to 0,$ $\displaystyle C^{\prime}_{G}(N):=\left\|\left\|G(t,x,y)-G_{N}(t,x,y)\right\|_{ L^{s^{\prime}}_{y}}\right\|_{L^{r}_{t}(0,T;L^{s}_{x})}\to 0$ $$

as $N\to\infty$,
respectively, where

$$ $G_{N}(t-\tau,x,y):=\sum_{m,n\in\Lambda_{N}}c_{n,m}\psi_{m}(\tau,y)\varphi_{n}( t,x),\quad 0\leq\tau<t<T,\ x,y\in D.$ (8) $$

###### Assumption 3 .

### 4.1 Main result

Our main result is the following quantitative approximation theorem.

###### Theorem 1 .

Suppose that $\mathcal{L}$ and $F$ satisfy
Assumptions 1 and 2, respectively, and the parameters $r,s$ satisfy the condition (3).
Then for any $R>0$, there exists a time $T>0$ such that the following statement holds under Assumption 3:
For any $\epsilon\in(0,1)$, there exist a depth $L$, the number of neurons $H$, a rank $N$, and $\Gamma\in\mathcal{NO}^{L,H,ReLU}_{N,\varphi,\psi}$ such that

$$ $\sup_{u_{0}\in B_{L^{\infty}}(R)}\|\Gamma^{+}(u_{0})-\Gamma(u_{0})\|_{L^{r}(0, T;L^{s})}\leq\epsilon.$ $$

Moreover, $L=L(\Gamma)$ and $H=H(\Gamma)$ satisfy

$$ $L(\Gamma)\leq C(\log(\epsilon^{-1}))^{2}\text{~{}~{}and~{}~{}}H(\Gamma)\leq C \epsilon^{-1}(\log(\epsilon^{-1}))^{2},$ $$

where $C>0$ is a constant depending on $\nu,M,F,D,p,d,r,s,R$, and $\mathcal{L}$, and $N=N(\Gamma)$ satisfies

$$ $C_{G}(N(\Gamma))\leq\epsilon\text{~{}~{}and~{}~{}}C^{\prime}_{G}(N(\Gamma)) \leq\epsilon.$ $$

Let us give a remark on model complexity of our neural operator $\Gamma$ in Theorem 1.
The model complexity depends on the depth $L(\Gamma)$, the number of neurons $H(\Gamma)$, and the rank $N$. As to $L(\Gamma)$ and $H(\Gamma)$,
we observe that they increase only squared logarithmically at most as $\epsilon\to 0$.
This is reasonable since each forward pass of the neural operator’s layers corresponds to a step in Picard’s iteration, and the convergence rate of Picard’s iteration decays exponentially (see Proposition 2 (ii)).
As to the rank $N$, Theorem 1 does not mention anything regarding the rates of the errors $C_{G}(N)$ and $C^{\prime}_{G}(N)$, but these error rates are also necessary for quantitatively estimating the model complexity of $\Gamma$. These error rates strongly depend on the choice of $\varphi$ and $\psi$, and the research results are enormous (see e.g. [, , , ]).
For example, in wavelet bases, the error rate is $O(N^{-s/d})$, where $s>0$ is the regularity of target function and $d$ is the spatial dimension (see [, Theorem 4.3.2]).
Therefore, if the Green function $G$ has a certain regularity, we may obtain a quantitative estimate for the rank $N$ such as the polynomial rate of $C_{G}(N)$ and $C^{\prime}_{G}(N)$, and as a result, it is shown that the model complexity of $\Gamma$ does not grow exponentially.
This prevents the exponential growth in model complexity, which often arises in general operator approximation.

###### Sketch of proof.

(See Appendix D for the exact proof).
The proof is based on the following scheme:

$$ $\displaystyle u_{0}$ $\displaystyle\xrightarrow[]{\ S_{\mathcal{L}}\ \,}u^{(1)}\xrightarrow[]{\ \ \, \Phi\ \ \ }u^{(2)}\xrightarrow[]{\ \ \,\Phi\ \ \ }\cdots\xrightarrow[]{\ \ \, \Phi\ \ \ }u^{(J)}(\cdots\xrightarrow[]{J\to\infty}u)$ (Picard’s iteration) $\displaystyle u_{0}$ $\displaystyle\xrightarrow[]{S_{\mathcal{L},N}}\hat{u}^{(1)}\xrightarrow[]{\Phi _{N,net}}\hat{u}^{(2)}\xrightarrow[]{\Phi_{N,net}}\cdots\xrightarrow[]{\Phi_{N ,net}}\hat{u}^{(J)}$ (Neural operator iteration) $$

where $S_{\mathcal{L}}$, $S_{\mathcal{L},N}$, $\Phi$, and $\Phi_{N,net}$ are defined as

$$ $\displaystyle S_{\mathcal{L}}[u_{0}](t,x):=\int_{D}G(t,x,y)u_{0}(y)\,dy,\quad S _{\mathcal{L},N}[u_{0}](t,x):=\int_{D}G_{N}(t,x,y)u_{0}(y)\,dy,$ $\displaystyle\Phi[u](t,x):=S_{\mathcal{L}}[u_{0}](t,x)+\int_{0}^{t}\int_{D}G(t ,x,y)F(u(\tau,y))\,dyd\tau,$ $\displaystyle\Phi_{N,net}[u](t,x):=S_{\mathcal{L},N}[u_{0}](t,x)+\int_{0}^{t} \int_{D}G_{N}(t,x,y)F_{net}(u(\tau,y))\,dyd\tau.$ $$

First, by the truncation of Picard’s iteration, the solution $u=\Gamma^{+}(u_{0})$ to (P’) is approximated by

$$ $u\approx u^{(J)}=\Phi^{[J-1]}\circ S_{\mathcal{L}}[u_{0}],$ $$

for sufficiently large $J$.
Next, by approximating the Green function $G$ and the nonlinearity $F$ with a truncated expansion $G_{N}$ and a neural network $F_{net}$, respectively, the contraction mapping $\Phi$ is approximated by $\Phi_{N,net}$ for sufficiently large $N$ (Lemmas 2 and  3), leading to

$$ $u^{(J)}\approx\hat{u}^{(J)}=\Phi^{[J-1]}_{N,net}\circ S_{\mathcal{L},N}[u_{0}],$ $$

(Lemma 4).
Finally, we represent $\hat{u}^{(J)}$ as a neural operator $\Gamma$ (Lemma 5).

The key idea is to mimic the Picard’s iteration.
The convergence rate of Picard’s iteration is given in Proposition 2 (ii), which reads $J\approx\log(\epsilon^{-1})$.
By applying the one dimensional case of [, Theorem 4.1], we approximate $F\approx F_{net}$ with the rates $L(F_{net})\lesssim(\log(\epsilon^{-1}))$ and $H(F_{net})\lesssim\epsilon^{-1}(\log(\epsilon^{-1}))$. Note that the curse of dimensionality that occurs in conventional neural network approximations does not appear here. Consequently, the rates of the depth $L(\Gamma)$ and the number of neurons $H(\Gamma)$ are evaluated as

$$ $\displaystyle L(\Gamma)\lesssim L(F_{net})\cdot J\approx\log\left(\epsilon^{-1 }\right)\cdot\log(\epsilon^{-1})=(\log(\epsilon^{-1}))^{2},$ $\displaystyle H(\Gamma)\lesssim H(F_{net})\cdot J\lesssim\epsilon^{-1}\left( \log(\epsilon^{-1})\right)\cdot\left(\log(\epsilon^{-1})\right)\lesssim \epsilon^{-1}\left(\log(\epsilon^{-1})\right)^{2},$ $$

where the first inequality is given by the exact construction of the neural operator $\Gamma$ discussed in the proof of Lemma 5 as

$$ $\displaystyle\Gamma(u_{0})=\tilde{W}^{\prime}\circ\left[\underbrace{(\tilde{W} +\tilde{K}_{N})\circ\tilde{F}_{net}\circ\cdots\circ(\tilde{W}+\tilde{K}_{N}) \circ\tilde{F}_{net}}_{J-1}\right]\circ(\tilde{K}^{(0)}_{N}+\tilde{b}^{(0)}_{N })(u_{0}),$ (9) $$

where $\tilde{W}^{\prime}$ and $\tilde{W}$ are some weight matrices, $\tilde{K}_{N}$ and $\tilde{K}^{(0)}_{N}$ are some non-local operators, $\tilde{b}^{(0)}_{N}$ is some bias function, and $\tilde{F}_{net}$ is some neural network.
∎

###### Remark 3 .

It is important to note that the hidden layers of the neural operator $\Gamma$ constructed in (9) are structured to repeatedly apply the same operator, $(\tilde{W}+\tilde{K}_{N})\circ\tilde{F}_{net}$, which significantly reduces memory consumption.
Our constructed neural operator is closely related to the deep equilibrium models [], which utilize weight-tied architectures to find fixed points.
More recently, inspired by deep equilibrium models, weight-tied neural operators have been experimentally investigated in the context of solving steady-state PDEs [].
Based on the above findings, our results also give weight-tied neural operators for (unsteady-state) parabolic PDEs and theoretically guarantee their approximation ability.

###### Remark 4 .

One of the next issues in approximation theory is to verify in what sense the obtained approximate solution is a good approximation.
In particular, it is important to discuss that an approximator $\Gamma$ obtained in Theorem 1 preserves desired properties in parabolic PDEs, such as the maximum principle and the comparison principle.
In the special case of the nonlinearity, we can show that approximator $\Gamma$ preserves the positivity. For more details, see Appendix F.

###### Theorem 1 .

###### Sketch of proof.

###### Remark 3 .

###### Remark 4 .

### 4.2 Approximation of long time solutions

Extension of a short time solution to a long time solution.
By Proposition 2 (i), for any initial data $u_{0}\in B_{L^{\infty}}(R)$, there exists a short time solution $u=\Gamma^{+}(u_{0})$ to (P’) on $[0,T]\times D$.
Then, if $u(T)$ at $t=T$ belongs to $B_{L^{\infty}}(R)$, the same solution operator $\Gamma^{+}$ can be applied to
$u(T)$ as an initial data, and hence, the solution can be extended to that on $[0,2T]\times D$.
Thus, as long as the supremum value of the solution is less than or equal to $R$, the solution can be extended with respect to $t$ by applying $\Gamma^{+}$ as many times as possible.
We denote by $u_{\max}$ the extended maximal solution with respect to $t$ and by $\kappa_{\max}=\kappa_{\max}(u_{0})\in\mathbb{N}\cup\{\infty\}$ the maximal number of extensions by $\Gamma^{+}$.
Then we can write $u_{\max}$ as

$$ $u_{\max}(t)={\Gamma^{+}}^{[\kappa]}(u_{0})(t)$ $$

for $(\kappa-1)T<t\leq\kappa T$ and $\kappa=1,\ldots,\kappa_{\max}$.

Approximation of the long time solution $u_{\max}$ by our neural operator.
Similarly, for our neural operator $\Gamma$ given in Theorem 1, the output function of $\Gamma$ can be extended with respect to time $t$ by applying $\Gamma$ as many times as possible.
We denote by $\hat{u}_{\max}$ the extended maximal output function of $\Gamma$ with respect to $t$ and by
$\hat{\kappa}_{\max}=\hat{\kappa}_{\max}(u_{0})\in\mathbb{N}\cup\{\infty\}$ the maximal number of extensions by $\Gamma$.
Then we can also write $\hat{u}_{\max}$ as

$$ $\hat{u}_{\max}(t):=\Gamma^{[\kappa]}(u_{0})(t)$ $$

for $(\kappa-1)T<t\leq\kappa T$ and $\kappa=1,\ldots,\hat{\kappa}_{\max}$.
Then the error between $u_{\max}$ and $\hat{u}_{\max}$ can be controlled by the error $\epsilon$ in Theorem 1 and another error $\tilde{\epsilon}$ given by

$$ $\tilde{\epsilon}:=\sum_{\kappa=1}^{\kappa^{*}-1}\|{\Gamma^{+}}^{[\kappa]}(u_{0 })(\kappa T)-\Gamma^{[\kappa]}(u_{0})(\kappa T)\|_{L^{\infty}}.$ $$

More precisely, we have the following result on the error between $u_{\max}$ and $\hat{u}_{\max}$.

###### Corollary 1 .

Let $\epsilon\in(0,1)$ and $R>0$. Suppose the same assumptions as in Theorem 1 and $\kappa^{*}:=\min\{\kappa_{\max},\hat{\kappa}_{\max}\}$ is a finite number.
Then there exists a constant $C>0$ such that

$$ $\|u_{\max}-\hat{u}_{\max}\|_{L^{r}(0,\kappa^{*}T;L^{s})}\leq C(\epsilon+\tilde {\epsilon})$ $$

for any $u_{0}\in B_{L^{\infty}}(R)$, where $T$ and $\Gamma$ are the same as in Theorem 1.

The proof is given in Appendix E.

###### Corollary 1 .

## 5 Discussion

A limitation of our method is that it applies only to parabolic PDEs that can be solved using the Banach’s fixed point theory. However, our approach is not restricted to parabolic PDEs and can be adapted to a broader class of nonlinear PDEs, because the Banach’s fixed point theory is widely applicable to many nonlinear PDEs. In addition, our proof based on Picard’s iteration is constructive, and this constructive approach may offer valuable insights for future experimental studies, particularly when incorporating constraints specific to the architecture of PDE tasks (for instance, weight-tied architecture discussed in Remark 3).
In future work, we will develop a theoretical framework applicable to a broader class of nonlinear PDEs, based on the ideas presented in this paper.

In numerical analysis, it is crucial that the solutions of discrete equations retain the properties of the original continuous solutions. Similarly, neural operators must also preserve these properties, especially when used as surrogate models of simulators in practical applications. While we have shown positivity preserving in a specific case of nonlinearity (Remark 4), future work may extend this result to more general nonlinear settings. Additionally, other important properties, such as the comparison principle and the smoothing effect, which are essential in parabolic equations, will be the focus of future studies.

#### Acknowledgments

T. Furuya is supported by JSPS KAKENHI Grant Number JP24K16949 and JST ASPIRE JPMJAP2329.
K. Taniguchi is supported by JSPS KAKENHI Grant Number 24K21316.

## Appendix

## Appendix A Examples of ℒ ℒ \mathcal{L} caligraphic_L (Assumption 1 )

This appendix provides some typical examples of $\mathcal{L}$ satisfying Assumption 1.

### A.1 Dirichlet Laplacian

Let $D$ be a bounded domain in $\mathbb{R}^{d}$ with $d\in\mathbb{N}$.
Let $\mathcal{L}=-\Delta_{\rm Dir}$ denote the Dirichlet Laplacian on $L^{2}(D)$.
Then the semigroup $\{S_{-\Delta_{\rm Dir}}(t)\}_{t\geq 0}$ corresponds to the solution operator of the initial boundary value problem of the linear heat equation

$$ $\begin{cases}\partial_{t}u-\Delta u=0\quad&\text{in }(0,T)\times D,\\ u(0)=u_{0}\quad&\text{in }D,\\ u=0\quad&\text{on }(0,T)\times\partial D,\end{cases}$ $$

where $\partial D$ denotes the boundary of $D$. The integral kernel $G_{\rm Dir}(t,x,y)$ of $S_{-\Delta_{\rm Dir}}(t)$ has the Gaussian upper bound

$$ $0\leq G_{\rm Dir}(t,x,y)\leq(4\pi t)^{-\frac{d}{2}}\exp\left(-\frac{|x-y|^{2}} {4t}\right),$ $$

for any $t>0$ and almost everywhere $x,y\in D$ (see e.g. [, Proposition 3.1]).
This bound yields
that $\mathcal{L}=-\Delta_{\rm Dir}$ satisfies Assumption 1 with $\nu=d/2$.

### A.2 Neumann or Robin Laplacian

Let $D$ be a bounded Lipschitz domain in $\mathbb{R}^{d}$ with $d\geq 1$. We consider the Robin Laplace $-\Delta_{\gamma}$ on $L^{2}(D)$ associated with a quadratic form

$$ $q_{\gamma}(f,g)=\int_{\Omega}\nabla f\cdot\nabla g\,dx+\int_{\partial\Omega} \gamma fg\,dS,$ $$

for any $f,g\in H^{1}(D)$, where $\gamma$ is a function $\partial\Omega\to\mathbb{R}$.
Note that $\gamma$ means a Robin boundary condition, and particularly, the case of $\gamma=0$ means the zero Neumann boundary condition. Hence, $-\Delta_{0}$ is the Neumann
Laplacian $-\Delta_{\rm Neu}$ on $L^{2}(D)$.
Assume that $\gamma\in L^{\infty}(\partial D)$ and $\gamma\geq 0$.
We denote by $G_{\gamma}(t,x,y)$ and $G_{\rm Neu}(t,x,y)$ the integral kernels of semigroups generated by $-\Delta_{\gamma}$ and $-\Delta_{\rm Neu}$, respectively.
Then $G_{\gamma}(t,x,y)$ and $G_{\rm Neu}(t,x,y)$ have the Gaussian upper bounds:

$$ $G_{\rm Dir}(t,x,y)\leq G_{\gamma}(t,x,y)\leq G_{\rm Neu}(t,x,y)\leq Ct^{-\frac {d}{2}}\exp\left(-\frac{|x-y|^{2}}{C^{\prime}t}\right),$ $$

for any $0<t\leq 1$ and almost everywhere $x,y\in D$ and for some $C,C^{\prime}>0$.
The first and second inequalities follow from domination of semigroups (see e.g. [, Theorem 2.24]),
and the proof of the last inequality can be found in [].
These bounds yield
that $-\Delta_{\rm Neu}$ and $-\Delta_{\gamma}$ satisfy Assumption 1 with $\nu=d/2$.

### A.3 Schrödinger operator with a singular potential

Let $D$ be a bounded domain in $\mathbb{R}^{d}$ with $d\geq 1$. We consider
the Schrödinger operator $\mathcal{L}=-\Delta+V$ with the zero Dirichlet boundary condition on $\partial D$, where
$V=V(x)$ is a real-valued measurable function on $D$ such that

$$ $V=V_{+}-V_{-},\quad V_{\pm}\geq 0,\quad V_{+}\in L^{1}_{\mathrm{loc}}(D)\quad \text{and}\quad V_{-}\in K_{d}(D).$ $$

We say that $V_{-}$ belongs to the Kato class $K_{d}(D)$ if

$$ $\displaystyle\left\{\begin{aligned} &\lim_{r\rightarrow 0}\sup_{x\in D}\int_{D \cap\{|x-y|<r\}}\frac{V_{-}(y)}{|x-y|^{d-2}}\,dy=0&\text{for }d\geq 3,\\ &\lim_{r\rightarrow 0}\sup_{x\in D}\int_{D\cap\{|x-y|<r\}}\log(|x-y|^{-1})V_{- }(y)\,dy=0&\text{for }d=2,\\ &\sup_{x\in D}\int_{D\cap\{|x-y|<1\}}V_{-}(y)\,dy<\infty&\text{for }d=1\,\end{ aligned}\right.,$ $$

(see [, Section A.2]).
It is readily seen that
the positive part $V_{+}$ is almost unrestricted, and the negative part $V_{-}$ allows for the singularity such as
$V_{-}(x)=1/|x|^{\alpha}$ with $0\leq\alpha<2$ if $d\geq 2$ and $0\leq\alpha<1$ if $d=1$.
Then it is known that
the integral kernel $G_{V}$ of semigroup generated by
the Schrödinger operator $\mathcal{L}=-\Delta+V$ satisfies
the Gaussian upper bound

$$ $|G_{V}(t,x,y)|\leq Ct^{-\frac{d}{2}}\exp\left(-\frac{|x-y|^{2}}{C^{\prime}t} \right),$ (10) $$

for any $0<t\leq 1$ and almost everywhere $x,y\in D$ and for some $C,C^{\prime}>0$ (see e.g. [, Proposition 3.1]).
Therefore, $\mathcal{L}=-\Delta+V$ satisfies
Assumption 1 with $\nu=d/2$.

### A.4 Uniformlly elliptic operator

We consider the following linear parabolic PDEs

$$ $\begin{cases}\partial_{t}u+\mathcal{L}u=0\quad&\text{in }(0,T)\times D,\\ u(0)=u_{0}\quad&\text{in }D,\end{cases}$ (11) $$

where

$$ $\mathcal{L}:=\sum_{k,j=1}^{d}\partial_{x_{k}}\left(a_{kj}\partial_{x_{j}} \right)+\sum_{k=1}^{d}b_{k}\partial_{x_{k}}+c_{0},$ $$

with real-valued functions $a_{k,j},b_{k},c_{0}$ on $D$ and
with domain $\mathcal{D}(\mathcal{L})$ that is a closed subspace of $H^{1}(D)$.
Note that the choice of domain $\mathcal{D}(\mathcal{L})$ corresponds to the choice of boundary condition.
For example, $\mathcal{D}(\mathcal{L})=H^{1}_{0}(D)$ and $H^{1}(D)$ mean the zero Dirichlet boundary condition and the zero Neumann boundary condition, respectively.
We assume that $\mathcal{L}$ satisfies the uniformly elliptic condition:

- •
$a_{k,j},b_{k},c_{0}\in L^{\infty}(D)\text{ for all }1\leq k,j\leq d$;
- •
There exists a constant $C>0$ such that
$\sum_{k,j=1}^{d}a_{k,j}(x)\xi_{k}\xi_{j}\geq C|\xi|^{2},$
for any $\xi\in\mathbb{R}^{d}$ and almost everywhere $x\in D$.

In addition, we also assume the following on $\mathcal{D}(\mathcal{L})$:

- •
$V$ is continuous embedded into $L^{2^{*}}(D)$, where $2^{*}=\frac{2d}{d-2}$ if $d\geq 3$, $2^{*}$ is any number in $(2,\infty)$ if $d=2$, and $2^{*}=\infty$ if $d=1$;
- •
If $u\in\mathcal{D}(\mathcal{L})$, then $\max\{0,u\},\min\{1,|u|\}\operatorname{sign}u\in\mathcal{D}(\mathcal{L})$;
- •
If $u\in\mathcal{D}(\mathcal{L})$, then $e^{\eta}u\in\mathcal{D}(\mathcal{L})$ for any real-valued function $\eta\in C^{\infty}(D)$ such that $\eta$ and $|\nabla\eta|$ are bounded on $D$.

Then it is known that the integral kernel of semigroup generated by $\mathcal{L}$ satisfies
the same Gaussian upper bound as (10) for any $0<t\leq 1$ and almost everywhere $x,y\in D$ (see [, Theorem 6.10]).
Therefore, the uniformly elliptic operator $\mathcal{L}$ satisfies
Assumption 1 with $\nu=d/2$.

### A.5 Other examples

Various other operators, such as the fractional Laplacian $(-\Delta)^{\alpha/2}$ ($\nu=d/\alpha$), the higher-order elliptic operator ($\nu=d/m$ with the highest order $m$) and the Schrödinger operator with a Dirac delta potential ($\nu=1/2$ with $d=1$), are possible, but these will not be mentioned here. For more information, see e.g. the books [, ] and the papers [, , , ]
and references therein.

## Appendix B Proof of Proposition 1

Let $D$ be a bounded domain in $\mathbb{R}^{d}$ with $d\in\mathbb{N}$.
We consider the Cauchy problem of nonlinear parabolic PDEs

$$ $\begin{cases}\partial_{t}u+\mathcal{L}u=F(u)\quad&\text{in }(0,T)\times D,\\ u(0)=u_{0}\quad&\text{in }D.\end{cases}$ (P) $$

This problem can be formally rewritten as the integral form

$$ $u(t)=S_{\mathcal{L}}(t)u_{0}+\int_{0}^{t}S_{\mathcal{L}}(t-\tau)F(u(\tau))\,d\tau.$ (P’) $$

Throughout this appendix, we always assume that $\mathcal{L}$ and $F$ satisfy Assumptions 1 and 2, respectively.
We study local well-posedness (LWP) for the integral equation (P’).
The LWP for (P’) means the existence of local in time solution, uniqueness of the solution, and continuous dependence on initial data for (P’), and
it has long been extensively studied (cf. [, , ]).
In this appendix, we provide the LWP results
that can be used as the basis of our study of neural operator and its approximation theorem.
We use the space $L^{r}(0,T;L^{s}(D))$ as a solution space.
For convenience, we set

$$ $\alpha=\alpha(\nu,p,s):=-\frac{\nu(p-1)}{s},\quad\beta=\beta(p,r):=\frac{r-p+1 }{r}.$ (12) $$

We have the following result.

###### Proposition 3 .

Let $d\in\mathbb{N}$ and $p\in(1,\infty)$ and $\nu>0$. Assume that $r,s$ satisfy

$$ $r,s\in[p,\infty]\quad\text{and}\quad\frac{\nu}{s}+\frac{1}{r}<\frac{1}{p-1}.$ (13) $$

Let $M>0$, $0<T\leq 1$ and $u_{0}\in L^{\infty}(D)$ be such that

$$ $\rho+\delta M\leq M,$ (14) $$

where $\rho$ and $\delta$ are defined by

$$ $\rho=\rho(T,\|u_{0}\|_{L^{\infty}}):=C_{\mathcal{L}}|D|^{\frac{1}{s}}T^{\frac{ 1}{r}}\|u_{0}\|_{L^{\infty}},$ $$

$$ $\delta=\delta(T,M):=2(\alpha\beta^{-1}+1)^{-\beta}C_{\mathcal{L}}C_{F}T^{ \alpha+\beta}M^{p-1},$ $$

respectively.
Define the map $\Phi=\Phi_{u_{0}}$ by

$$ $\Phi[u](t):=S_{\mathcal{L}}(t)u_{0}+\int_{0}^{t}S_{\mathcal{L}}(t-\tau)F(u( \tau))\,d\tau,$ $$

for $t\in[0,T]$ and the complete metric space $X=X(T,M)$ by

$$ $X:=B_{L^{r}(0,T;L^{s})}(M)=\{u\in L^{r}(0,T;L^{s}(D)):\|u\|_{L^{r}(0,T;L^{s})} \leq M\},$ $$

equipped with the metric

$$ $\mathrm{d}(u,v):=\|u-v\|_{L^{r}(0,T;L^{s})}.$ $$

Then $\Phi$ is a contraction mapping from $X$ to itself.
Consequently, there exists a unique function $u\in X$ such that $\Phi[u]=u$.

###### Proof.

The proof is based on the standard fixed point argument.
Let us take the parameters $r,s$ satisfying

$$ $r,s\in[1,\infty],\quad\frac{s}{p}\geq 1,\quad\frac{r}{p}\geq 1,\quad 0<\beta \leq 1\quad\text{and}\quad\alpha+\beta>0.$ (15) $$

Take $M>0$, $0<T\leq 1$ and $u_{0}\in L^{\infty}(D)$ such that

$$ $\rho+\delta M\leq M\quad\text{and}\quad\delta<1.$ (16) $$

Let $u\in X$. Then

$$ $\begin{split}\|\Phi[u]\|_{L^{r}(0,T;L^{s})}&\leq\|S_{\mathcal{L}}(t)u_{0}\|_{L ^{r}(0,T;L^{s})}+C_{F}\left\|\int_{0}^{t}S_{\mathcal{L}}(t-\tau)(|u(\tau)|^{p- 1}u(\tau))\,d\tau\right\|_{L^{r}(0,T;L^{s})}\\ &=:I+II.\end{split}$ $$

As for the first term,
it follows from Hölder’s inequality and (1) that

$$ $I\leq|D|^{\frac{1}{s}}T^{\frac{1}{r}}\|S_{\mathcal{L}}(t)u_{0}\|_{L^{\infty}(0 ,T;L^{\infty})}\leq C_{\mathcal{L}}|D|^{\frac{1}{s}}T^{\frac{1}{r}}\|u_{0}\|_{ L^{\infty}}=\rho.$ $$

As for the second term, it follows from (1) that

$$ $\begin{split}II&\leq C_{\mathcal{L}}C_{F}\left\|\int_{0}^{t}(t-\tau)^{-\nu( \frac{p}{s}-\frac{1}{s})}\||u(\tau)|^{p-1}u(\tau)\|_{L^{\frac{s}{p}}}\,d\tau \right\|_{L^{r}((0,T))}\\ &\leq C_{\mathcal{L}}C_{F}\left\|\int_{0}^{t}(t-\tau)^{\alpha}\|u(\tau)\|_{L^{ s}}^{p}\,d\tau\right\|_{L^{r}((0,T))}\\ &\leq C_{\mathcal{L}}C_{F}\|f*g\|_{L^{r}(\mathbb{R})},\end{split}$ $$

where $s\geq 1$ and $s/p\geq 1$ are required and the functions $f,g$ are defined by

$$ $f(t)=f_{T}(t):=\begin{cases}t^{\alpha}\quad&\text{for }0<t\leq T,\\ 0&\text{otherwise},\end{cases}$ $$

and

$$ $g(t)=g_{T}(t):=\begin{cases}\|u(t)\|_{L^{s}}^{p}\quad&\text{for }0<t\leq T,\\ 0&\text{otherwise}.\end{cases}$ $$

By Young’s inequality, we have

$$ $\begin{split}\|f*g\|_{L^{r}(\mathbb{R})}&\leq\|f\|_{L^{\frac{1}{\beta}}( \mathbb{R})}\|g\|_{L^{\frac{r}{p}}(\mathbb{R})}\\ &=\left(\alpha\beta^{-1}+1\right)^{-\beta}T^{\alpha+\beta}\|u\|_{L^{r}(0,T;L^{ s})}^{p}\\ &\leq\left(\alpha\beta^{-1}+1\right)^{-\beta}T^{\alpha+\beta}M^{p},\\ \end{split}$ $$

where we require $r/p\geq 1$, $0<\beta\leq 1$ and
$\alpha+\beta>0$ and we calculate

$$ $\|f\|_{L^{\frac{1}{\beta}}(\mathbb{R})}=\left(\alpha\beta^{-1}+1\right)^{- \beta}T^{\alpha+\beta}.$ $$

Summarizing the estimates obtained above, and then using (14), we have

$$ $\|\Phi[u]\|_{L^{r}(0,T;L^{s})}\leq\rho+C_{\mathcal{L}}C_{F}\left(\alpha\beta^{ -1}+1\right)^{-\beta}T^{\alpha+\beta}M^{p}=\rho+\delta M\leq M,$ $$

which implies that $\Phi$ is a mapping from $X$ to itself.

Next, we estimate the metric $\mathrm{d}(\Phi[u_{1}],\Phi[u_{2}])$.
Let $u_{1},u_{2}\in X$.
By using (2), (1) and Hölder’s inequality, we have

$$ $\begin{split}\mathrm{d}(\Phi[u_{1}],\Phi[u_{2}])&\leq C_{\mathcal{L}}C_{F} \left\|\int_{0}^{t}(t-\tau)^{\alpha}\|\max_{i=1,2}|u_{i}(\tau)|^{p-1}|u_{1}( \tau)-u_{2}(\tau)|\|_{L^{\frac{s}{p}}}\,d\tau\right\|_{L^{r}((0,T))}\\ &\leq 2C_{\mathcal{L}}C_{F}\max_{i=1,2}\left\|\int_{0}^{t}(t-\tau)^{\alpha}\|u _{i}(\tau)\|_{L^{s}}^{p-1}\|u_{1}(\tau)-u_{2}(\tau)\|_{L^{s}}\,d\tau\right\|_{ L^{r}((0,T))}\\ &\leq 2C_{\mathcal{L}}C_{F}\max_{i=1,2}\|f*h_{i}\|_{L^{r}(\mathbb{R})},\end{split}$ $$

where $f$ is the same function as above and $h_{i}$ is defined by

$$ $h_{i}(t)=h_{i,T}(t):=\begin{cases}\|u_{i}(\tau)\|_{L^{s}}^{p-1}\|u_{1}(\tau)-u _{2}(\tau)\|_{L^{s}}\quad&\text{for }0<t\leq T,\\ 0&\text{otherwise}.\end{cases}$ $$

By Young’s inequality and Hölder’s inequality, we have

$$ $\begin{split}\mathrm{d}(\Phi[u_{1}],\Phi[u_{2}])&\leq 2C_{\mathcal{L}}C_{F} \max_{i=1,2}\|f*h_{i}\|_{L^{r}(\mathbb{R})}\\ &\leq 2C_{\mathcal{L}}C_{F}\max_{i=1,2}\|f\|_{L^{\frac{1}{\beta}}(\mathbb{R})} \|h_{i}\|_{L^{\frac{r}{p}}(\mathbb{R})}\\ &\leq 2C_{\mathcal{L}}C_{F}\left(\alpha\beta^{-1}+1\right)^{-\beta}T^{\alpha+ \beta}\max_{i=1,2}\|u_{i}\|_{L^{r}(0,T;L^{s})}^{p-1}\|u_{1}-u_{2}\|_{L^{r}(0,T ;L^{s})}\\ &\leq 2C_{\mathcal{L}}C_{F}\left(\alpha\beta^{-1}+1\right)^{-\beta}T^{\alpha+ \beta}M^{p-1}\mathrm{d}(u_{1},u_{2})=\delta\mathrm{d}(u_{1},u_{2}).\end{split}$ $$

From the second inequality $\delta<1$ in (16), we show that
$\Phi$ is a contraction mapping from $X$ into itself.
Therefore, Banach’s fixed point theorem allows us to prove that there exists uniquely
a function $u\in X$ such that $u=\Phi[u]$.
Finally, noting the condition (15) is equivalent to (13) and the condition (16) is also equivalent to (14), we conclude Proposition 3.
∎

Since $\min\{1/r,\alpha+\beta\}$ is always positive
under the condition (13),
for any $M>0$ and $u_{0}\in L^{\infty}(D)$, we can take $T>0$ sufficiently small so that (14) holds.
As a result, we have the following LWP result.

###### Corollary 2 (Proposition 1 ) .

Assume that $r,s$ satisfy
(13). Then,
for any $u_{0}\in L^{\infty}(D)$, there exist a time $T=T(u_{0})>0$ and a unique solution $u\in L^{r}(0,T;L^{s}(D))$ to (P’). Moreover, for any $u_{0},v_{0}\in L^{\infty}(D)$, the solutions $u,v$ to (P’) with $u(0)=u_{0}$ and $v(0)=v_{0}$ satisfy
the inequality

$$ $\|u-v\|_{L^{r}(0,T^{\prime};L^{s})}\leq C\|u_{0}-v_{0}\|_{L^{\infty}},$ $$

where $T^{\prime}=\min\{T(u_{0}),T(v_{0})\}$.

###### Proof.

The former immediately follows from Proposition 3. The latter is also shown in a similar way to the proof of Proposition 3. In fact, let $u,v\in X$ be solutions to (P’) with initial data $u_{0}$ and $v_{0}$, respectively. Then

$$ $\|u-v\|_{L^{r}(0,T;L^{s})}\leq C_{\mathcal{L}}|D|^{\frac{1}{s}}T^{\frac{1}{r}} \|u_{0}-v_{0}\|_{L^{\infty}}+\delta\|u-v\|_{L^{r}(0,T;L^{s})}.$ $$

If $T$ is sufficiently small so that $\delta<1$, then we obtain

$$ $\|u-v\|_{L^{r}(0,T;L^{s})}\leq\frac{C_{\mathcal{L}}|D|^{\frac{1}{s}}T^{\frac{1 }{r}}}{1-\delta}\|u_{0}-v_{0}\|_{L^{\infty}}.$ $$

The proof is finished.
∎

###### Proposition 3 .

###### Proof.

###### Corollary 2 (Proposition 1 ) .

###### Proof.

## Appendix C Applications to FNOs and WNOs

In this appendix, we apply Theorem 1 to FNOs and WNOs for some operator $\mathcal{L}$.
Let us recall that, to ensure our quantitative approximation theorem (Theorem 1), the families
$\varphi:=\{\varphi_{n}\}_{n}$ and $\psi:=\{\psi_{m}\}_{m}$ of functions need to approximate
the Green function $G$ of the linear equation $\partial_{t}u+\mathcal{L}u=0$ (i.e. the integral kernel $G$ of semigroup generated by $\mathcal{L}$) such that

$$ $G(t-\tau,x,y)=\sum_{m,n\in\Lambda}c_{m,n}\psi_{m}(\tau,y)\varphi_{n}(t,x),$ (17) $$

$$ $G_{N}(t-\tau,x,y)=\sum_{m,n\in\Lambda_{N}}c_{m,n}\psi_{m}(\tau,y)\varphi_{n}(t ,x)$ $$

for $0\leq\tau<t<T$ and $x,y\in D$,
where $\Lambda$ is an index set that is either finite or countably infinite and $\Lambda_{N}$ is a subset of $\Lambda$ with its cardinality $|\Lambda_{N}|=N\in\mathbb{N}$. The convergence (17) means the sense that

$$ $C_{G}(N):=\left\|\left\|G(t-\tau,x,y)-G_{N}(t-\tau,x,y)\right\|_{L^{r^{\prime} }_{\tau}(0,t;L^{s^{\prime}}_{y})}\right\|_{L^{r}_{t}(0,T;L^{s}_{x})}\to 0$ (18) $$

and

$$ $C^{\prime}_{G}(N):=\left\|\left\|G(t,x,y)-G_{N}(t,x,y)\right\|_{L^{s^{\prime}} _{y}}\right\|_{L^{r}_{t}(0,T;L^{s}_{x})}\to 0$ (19) $$

as $N\to\infty$ for some $r,s\in[1,\infty]$ satisfying the condition (3). The following conditions are required for the arguments in this appendix:

$$ $\left\|\left\|G(t-\tau,x,y)\right\|_{L^{r^{\prime}}_{\tau}(0,t;L^{s^{\prime}}_ {y})}\right\|_{L^{r}_{t}(0,T;L^{s}_{x})}<\infty,$ (20) $$

$$ $\left\|\left\|G(t,x,y)\right\|_{L^{s^{\prime}}_{y}}\right\|_{L^{r}_{t}(0,T;L^{ s}_{x})}<\infty.$ (21) $$

###### Remark 5 .

Let us give some remarks on (20) and (21).

- (a)
It is seen that
(20) and (21) always hold with $r=s=\infty$ under Assumption 1 (see Lemma 1 in Appendix D).
- (b)
For example, the operators $\mathcal{L}$ appeared in Appendix A.1–A.4 satisfy the Gaussian upper bound (10), and hence, the corresponding Green functions satisfy
(20) and (21) if
$\frac{d}{s}+\frac{2}{r}<2\quad\text{and}\quad\frac{d}{s}<\frac{2}{r}.$
(22)

###### Remark 5 .

### C.1 Fourier neural operators

Before stating FNOs, we provide a remark on partial sums of multiple Fourier series. From the perspective of convergence issues, the rectangular partial sum is preferable to the spherical partial sum. In fact,
for $f\in L^{1}(\mathbb{T}^{d})$, the Fourier coefficient of $f$ is defined by

$$ $c_{n}=\hat{f}(n):=\int_{\mathbb{T}^{d}}f(x)e^{-2\pi inx}\,dx,\quad n=(n_{1},n_ {2},\cdots,n_{d})\in\mathbb{Z}^{d},$ $$

where $\mathbb{T}^{d}=\mathbb{R}^{d}/\mathbb{Z}^{d}$
is the $d$-dimensional torus, $x=(x_{1},x_{2},\ldots,x_{d})$ and $nx:=n_{1}x_{1}+\cdots+n_{d}x_{d}$.
We define the rectangular partial sum of multiple Fourier series of $f$ by

$$ $S_{N}(f)(x):=\sum_{|n_{1}|,\cdots,|n_{d}|\leq N}c_{n}e^{2\pi inx},\quad N\in \mathbb{N}.$ $$

Then we have the following:

###### Theorem 2 ( [ 44 , 17 ] ) .

Let $q>1$. Then

$$ $\lim_{N\to\infty}S_{N}(f)(x)=f(x)\quad\text{a.e. $x$}\quad\text{and}\quad\lim_ {N\to\infty}S_{N}(f)=f\quad\text{in }L^{q}(\mathbb{T}^{d}),$ $$

for any $f\in L^{q}(\mathbb{T}^{d})$.

This theorem can be easily extended to periodic functions $f$ on a general $d$-dimensional rectangle with usual modifications. In contrast,
the spherical partial sum of multiple Fourier series of $f\in L^{q}(\mathbb{T}^{d})$ converges to $f$ in $L^{q}$ only if $q=2$, and for its pointwise convergence, it remains an open problem even when $q=2$.

Next, we mention the approximation theorem on FNOs.
Here, we apply the Fourier expansion with respect to both time $t$ and space $x$ to define FNOs. This type of FNO is used in the implementation to approximate some PDEs (see e.g. [, Section 7.3]).

For simplicity, we consider only the case $d=1$, $r=s=2$ and $\mathcal{L}$ is one of Appendixes A.1–A.4 (which implies $\nu=d/2$).
Then the condition (22) is satisfied and Theorem 1 can be applied to FNOs under the condition (3), i.e., $1<p<1+\frac{4}{d+2}$,
if we can show that (18) and (19) hold for the Fourier basis.

In fact, we consider
the nonlinear parabolic PDEs (P),
where $D$ is a bounded interval in $\mathbb{R}$. As $G$ is not a periodic function, a zero extension of $G$ is considered here, so that the multiple Fourier series can be applied.
We take $\tilde{D}$ as a bounded interval which includes $D$ and with its length $L$, and we
define
the zero extension $\tilde{G}$ of $G=\tilde{G}(t-\tau,x,y)$ to $[-T,2T]\times[-T,2T]\times\tilde{D}\times\tilde{D}$ by

$$ $\tilde{G}(t-\tau,x,y)=\begin{cases}G(t-\tau,x,y)\quad&\text{if }0\leq\tau<t<T, \ x,y\in D\\ 0&\text{otherwise}.\end{cases}$ $$

Then, thanks to (20) and (21), $\tilde{G}$ has the multiple Fourier series

$$ $\tilde{G}(t-\tau,x,y)=\lim_{N\to\infty}\sum_{m,n\in\Lambda_{N}}c_{m,n}e^{\frac {2\pi in_{1}t}{3T}}e^{\frac{2\pi im_{1}\tau}{3T}}e^{\frac{2\pi in_{2}x}{L}}e^{ \frac{2\pi im_{2}}{y}}$ $$

for almost everywhere in $(t,\tau,x,y)\in[-T,2T]\times[-T,2T]\times\tilde{D}\times\tilde{D}$ and it also converges
in $L^{2}([-T,2T]\times[-T,2T]\times\tilde{D}\times\tilde{D})$,
where $\Lambda_{N}:=\{(m,n)=(m_{1},m_{2},n_{1},n_{2})\in\mathbb{Z}^{4}:|m_{1}|,|m_{2}
|,|n_{1}|,|n_{2}|\leq N^{1/4}\}$ and $c_{n}$ is the Fourier coefficient given by

$$ $c_{m,n}:=\int_{-T}^{2T}\int_{-T}^{2T}\int_{\tilde{D}}\int_{\tilde{D}}\tilde{G} (t-\tau,x,y)e^{-\frac{2\pi in_{1}t}{3T}}e^{-\frac{2\pi in_{2}\tau}{3T}}e^{- \frac{2\pi in_{3}x}{L}}e^{-\frac{2\pi in_{4}y}{L}}\,dtd\tau dxdy.$ $$

Hence, putting
$\varphi_{n}(t,x):=e^{\frac{2\pi in_{1}t}{3T}}e^{\frac{2\pi in_{2}x}{L}}$ and
$\psi_{m}(\tau,y):=e^{\frac{2\pi im_{1}\tau}{3T}}e^{\frac{2\pi im_{2}y}{L}}$, we have the expansion (17) of $G$ with (18) and (19).
Therefore, by choosing these $\varphi=\{\varphi_{n}\}_{n}$ and $\psi=\{\psi_{m}\}_{m}$ in Definition 1, we have the following result on the FNO as a corollary of Theorem 1.

###### Corollary 3 .

Let $d=1$, $r=s=2$, $\mathcal{L}$ be one of Appendixes A.1–A.4 ($\nu=1/2$), $\varphi=\{\varphi_{n}\}_{n}$ and $\psi=\{\psi_{m}\}_{m}$ as above, and
$F$ satisfy Assumption 2 with $1<p<7/3$.
Then for any $R>0$, there exists a time $T>0$ such that the following statement holds:
For any $\epsilon\in(0,1)$, there exist a depth $L$, the number of neurons $H$, a rank $N$ and an FNO $\Gamma\in\mathcal{NO}^{L,H,ReLU}_{N,\varphi,\psi}$ such that

$$ $\sup_{u_{0}\in B_{L^{\infty}}(R)}\|\Gamma^{+}(u_{0})-\Gamma(u_{0})\|_{L^{2}(0, T;L^{2})}\leq\epsilon.$ $$

Moreover, $L=L(\Gamma)$ and $H=H(\Gamma)$ satisfy

$$ $L(\Gamma)\leq C(\log(\epsilon^{-1}))^{2}\text{~{}~{}and~{}~{}}H(\Gamma)\leq C \epsilon^{-1}(\log(\epsilon^{-1}))^{2},$ $$

where $C>0$ is a constant depending on $\nu,M,F,D,p,d,r,s,R$ and $\mathcal{L}$.

###### Theorem 2 ( [ 44 , 17 ] ) .

###### Corollary 3 .

### C.2 Haar wavelet neural operators

A wavelet is an oscillatory function with zero mean, used to analyze signals or functions locally at different scales.
There are various types of wavelets, such as Haar, Meyer, Morlet, and Shannon and the systems constructed by these wavelets are unconditional bases for $L^{q}(\mathbb{R}^{d})$ with $1<q<\infty$
(see e.g. [, ]). Wavelets can also be defined on domains other than on $\mathbb{R}^{d}$ (see e.g. []).
Several results on wavelets for mixed Lebesgue spaces are also known (see e.g. [, , ]).

In this appendix, we use the fact that the multiparameter Haar system is an unconditional basis for the mixed Lebesgue spaces $L^{\vec{q}}((0,1)^{d})$ with $\vec{q}=(q_{1},\ldots,q_{d})\in(1,\infty)^{d}$ by [, Theorem 4.9].

The set of dyadic intervals in $[0,1]$ is defined by

$$ $\mathcal{D}:=\left\{I_{j,k}:=\left[\frac{k}{2^{j-1}},\frac{k+1}{2^{j-1}}\right ):0\leq k<2^{j-1},\ j\in\mathbb{N}\right\}.$ $$

For $I\in\mathcal{D}$, we denote by $I_{\mathrm{left}}$ and $I_{\mathrm{right}}$
the left and right halves of $I$, respectively, and we define the Haar function $h_{I}$ by $h_{I}:=\chi_{I_{\mathrm{left}}}-\chi_{I_{\mathrm{right}}}$.
Here, $\chi_{I}$ is the characteristic function of an interval $I$.
The set of dyadic hyper-rectangles in $[0,1]^{d}$ is defined by

$$ $\mathcal{R}_{d}:=\left\{\mathcal{I}:=I_{1}\times I_{2}\times\cdots\times I_{d} :I_{1},I_{2},\cdots,I_{d}\in\mathcal{D}\right\}.$ $$

For $\mathcal{I}\in\mathcal{R}_{d}$, we define

$$ $h_{\mathcal{I}}(x):=h_{I_{1}}(x_{1})h_{I_{2}}(x_{2})\cdots h_{I_{d}}(x_{d})$ $$

for $x=(x_{1},x_{2},\ldots,x_{d})\in[0,1]^{d}$.
Then the multiparameter Haar system defined by $\{h_{\mathcal{I}}:\mathcal{I}\in\mathcal{R}_{d}\}$
is an unconditional basis for the mixed Lebesgue spaces $L^{\vec{q}}((0,1)^{d})$ with $\vec{q}=(q_{1},\ldots,q_{d})\in(1,\infty)^{d}$
(see [, Theorem 4.8]).

Let $\tilde{G}$ be the zero extension of the Green function $G$ to a function on $[0,1]^{2d+2}$ ($=[0,1]^{1+d+1+d}$) with respect to $t,x,\tau,y$.
Then
we can apply [, Theorem 4.8] to $\tilde{G}$ and then

$$ $\tilde{G}(t-\tau,x,y)=\lim_{N\to\infty}\sum_{\mathcal{I}\in\mathcal{R}_{2d+2,N }}c_{\mathcal{I}}h_{\mathcal{I}_{t}}(t)h_{\mathcal{I}_{x}}(x)h_{\mathcal{I}_{ \tau}}(\tau)h_{\mathcal{I}_{y}}(y)$ $$

in the sense of mixed Lebesgue norms of $L^{\vec{q}}((0,1)^{2d+2})$ with all $\vec{q}=(q_{1},\ldots,q_{d})\in(1,\infty)^{2d+2}$,
where $\mathcal{R}_{2d+2,N}$ is a subset of $\mathcal{R}_{2d+2}$ with its cardinality $|\mathcal{R}_{2d+2,N}|=N$,
$c_{\mathcal{I}}$ are wavelet coefficients, and
$\mathcal{I}=\mathcal{I}_{t}\times\mathcal{I}_{x}\times\mathcal{I}_{\tau}\times
\mathcal{I}_{y}\in\mathcal{D}\times\mathcal{D}^{d}\times\mathcal{D}\times
\mathcal{D}^{d}$.
Hence, putting
$c_{m,n}=c_{\mathcal{I}}$,
$\varphi_{n}(t,x)=h_{\mathcal{I}_{t}}(t)h_{\mathcal{I}_{x}}(x)$ and $\psi_{m}(\tau,y)=h_{\mathcal{I}_{\tau}}(\tau)h_{\mathcal{I}_{y}}(y)$ with $m=(k_{1},j_{1})\in(\mathbb{N}\cup\{0\})^{d+1}\times\mathbb{N}^{d}$ and $n=(k_{2},j_{2})\in(\mathbb{N}\cup\{0\})^{d+1}\times\mathbb{N}^{d}$, we have the expansion (17) of $G$ with (18) and (19), where $\Lambda_{N}$ is an index set of $(k_{1},j_{1},k_{2},j_{2})$ corresponding to $\mathcal{R}_{2d+2,N}$ with its cardinality $|\Lambda_{N}|=N$.
Therefore, by choosing these $\varphi=\{\varphi_{n}\}_{n}$ and $\psi=\{\psi_{m}\}_{m}$ in Definition 1, we have the following result on the Haar WNO as a corollary of Theorem 1.

###### Corollary 4 .

Let $d\in\mathbb{N}$. Suppose $\mathcal{L}$ is one of Appendixes A.1–A.4 ($\nu=d/2$) and $F$ satisfies Assumption 2.
Let $r,s\in(1,\infty)$ satisfy the conditions (3) and (22), and let
$\varphi=\{\varphi_{n}\}_{n}$ and $\psi=\{\psi_{m}\}_{m}$ be as above.
Then, for any $R>0$, there exists a time $T>0$ such that the following statement holds:
For any $\epsilon\in(0,1)$, there exist a depth $L$, the number of neurons $H$, a rank $N$ and a Haar WNO $\Gamma\in\mathcal{NO}^{L,H,ReLU}_{N,\varphi,\psi}$ such that

$$ $\sup_{u_{0}\in B_{L^{\infty}}(R)}\|\Gamma^{+}(u_{0})-\Gamma(u_{0})\|_{L^{r}(0, T;L^{s})}\leq\epsilon.$ $$

Moreover, $L=L(\Gamma)$ and $H=H(\Gamma)$ satisfy

$$ $L(\Gamma)\leq C(\log(\epsilon^{-1}))^{2}\text{~{}~{}and~{}~{}}H(\Gamma)\leq C \epsilon^{-1}(\log(\epsilon^{-1}))^{2},$ $$

where $C>0$ is a constant depending on $\nu,M,F,D,p,d,r,s,R$ and $\mathcal{L}$.

###### Corollary 4 .

## Appendix D Proof of Theorem 1

Suppose that $\mathcal{L}$ and $F$ satisfy
Assumptions 1 and 2, respectively, and the parameters $r,s$ satisfy the condition (3).
The conditions for the parameters $R,T,M,M^{\prime}$ that appear in the proof are specified here.
Let us take $R,M,M^{\prime}>0$ and $T\in(0,1]$ such that

$$ $\begin{dcases}C_{\mathcal{L}}|D|^{\frac{1}{s}}T^{\frac{1}{r}}R+(\alpha\beta^{- 1}+1)^{-\beta}C_{\mathcal{L}}C_{F}T^{\alpha+\beta}M^{p}\leq M,\\ 2C_{\mathcal{L}}R+2C_{\mathcal{L}}T(1+C_{F}M^{\prime p})\leq M^{\prime},\\ T^{\frac{1}{r}}|D|^{\frac{1}{s}}M^{\prime}\leq M,\end{dcases}$ (23) $$

where $\alpha$ and $\beta$ are the constants defined in (12).
The first condition in (23) ensures the LWP for (P’) (Propositions 1 and 2). See Proposition 3 and (14) in Appendix B for more details.
The second condition ensures that the approximate mappings $\Phi_{N_{k}}$ and $\Phi_{N_{k},net}$ in the proof map from $B_{L^{\infty}(0,T;L^{\infty})}(M^{\prime})$ into itself.
The third condition ensures $B_{L^{\infty}(0,T;L^{\infty})}(M^{\prime})\subset B_{L^{r}(0,T;L^{s})}(M)$.
We note to take the parameters in the following order: For any $R,M>0$, we can take a sufficiently large $M^{\prime}$ and a sufficiently small $T$ so that (23) holds.

To begin with, we prepare the following lemma.

###### Lemma 1 .

Under Assumption 1, we have

$$ $\displaystyle\left\|\left\|G(t-\tau,x,y)\right\|_{L^{1}_{\tau}(0,t;L^{1}_{y})} \right\|_{L^{\infty}_{t}(0,T;L^{\infty}_{x})}\leq C_{\mathcal{L}}T,$ $\displaystyle\left\|\left\|G(t,x,y)\right\|_{L^{1}_{y}}\right\|_{L^{\infty}_{t }(0,T;L^{\infty}_{x})}\leq C_{\mathcal{L}},$ $$

for any $T\in(0,1]$.

###### Proof.

From Assumption 1 and [, Lemma B.1], we see that

$$ $\displaystyle\left\|G(t-\tau,x,y)\right\|_{L^{\infty}_{x}L^{1}_{y}}=\|S_{ \mathcal{L}}(t-\tau)\|_{L^{\infty}\to L^{\infty}}\leq C_{\mathcal{L}},$ $$

which implies that

$$ $\begin{split}\left\|\left\|G(t-\tau,x,y)\right\|_{L^{1}_{\tau}(0,t;L^{1}_{y})} \right\|_{L^{\infty}_{t}(0,T;L^{\infty}_{x})}\leq C_{\mathcal{L}}T.\end{split}$ $$

Similarly, as

$$ $\displaystyle\left\|G(t,x,y)\right\|_{L^{\infty}_{x}L^{1}_{y}}=\|S_{\mathcal{L }}(t)\|_{L^{\infty}\to L^{\infty}}\leq C_{\mathcal{L}},$ $$

we also have

$$ $\left\|\left\|G(t,x,y)\right\|_{L^{1}_{y}(D)}\right\|_{L^{\infty}_{t}(0,T;L^{ \infty}_{x})}\leq C_{\mathcal{L}}.$ $$

The proof of Lemma 1 is finished.
∎

We divide the proof of Theorem 1 into four steps.

Step 1.
We define $\Phi_{N}$ by

$$ $\Phi_{N}[u](t,x):=\int_{D}G_{N}(t,x,y)u_{0}(y)dy+\int_{0}^{t}\int_{D}G_{N}(t- \tau,x,y)F(u(\tau,y))d\tau dy.$ $$

Since $C_{G}(N),C^{\prime}_{G}(N)\to 0$ as $N\to\infty$ with $r^{\prime},s^{\prime}\geq 1$ and $D\subset\mathbb{R}^{d}$ is bounded, we see that

$$ $\displaystyle\left\|\left\|G(t-\tau,x,y)-G_{N}(t-\tau,x,y)\right\|_{L^{1}_{ \tau}(0,t;L^{1}_{y})}\right\|_{L^{r}_{t}(0,T;L^{s}_{x})}\to 0,$ $\displaystyle\left\|\left\|G(t,x,y)-G_{N}(t,x,y)\right\|_{L^{1}_{y}}\right\|_{ L^{r}_{t}(0,T;L^{s}_{x})}\to 0,$ $$

as $N\to\infty$, and hence,
there exists a subsequence
$\{N_{k}\}_{k\in\mathbb{N}}\subset\mathbb{N}$ such that

$$ $\displaystyle\left\|G_{N_{k}}(t-\tau,x,y)\right\|_{L^{1}_{\tau}(0,t;L^{1}_{y}) }\to\left\|G(t-\tau,x,y)\right\|_{L^{1}_{\tau}(0,t;L^{1}_{y})},$ $\displaystyle\left\|G_{N_{k}}(t,x,y)\right\|_{L^{1}_{y}}\to\left\|G(t,x,y) \right\|_{L^{1}_{y}},$ $$

for almost everywhere $t\in(0,T)$ and $x\in D$ as
$k\to\infty$.
Then, for any $\epsilon>0$,
there exists $k_{\epsilon}\in\mathbb{N}$ such that for any $k\geq k_{\epsilon}$,

$$ $\displaystyle\left\|\left\|G_{N_{k}}(t-\tau,x,y)\right\|_{L^{1}_{\tau}(0,t;L^{ 1}_{y})}\right\|_{L^{\infty}_{t}(0,T;L^{\infty}_{x})}\leq 2C_{\mathcal{L}}T,$ $\displaystyle\left\|\left\|G_{N_{k}}(t,x,y)\right\|_{L^{1}_{y}}\right\|_{L^{ \infty}_{t}(0,T;L^{\infty}_{x})}\leq 2C_{\mathcal{L}},$ $\displaystyle C_{G}(N_{k}),C^{\prime}_{G}(N_{k})\leq\epsilon,$ $$

where the first and second inequalities follow from Lemma 1 and
the last one follows from Assumption 3.

###### Lemma 2 .

For any $u_{0}\in B_{L^{\infty}}(R)$ and $k\geq k_{\epsilon}$, the following statements hold:

- (i)
$\Phi_{N_{k}}$ is a map from $B_{L^{\infty}(0,T;L^{\infty})}(M^{\prime})$ to itself.
- (ii)
There exists a constant $C_{1}>0$ such that
$\|\Phi[u]-\Phi_{N_{k}}[u]\|_{L^{r}(0,T;L^{s})}\leq C_{1}\epsilon,$
for $u\in B_{L^{\infty}(0,T;L^{\infty})}(M^{\prime})$.
Here, $C_{1}=|D|^{\frac{1}{s}}R+C_{F}M^{p}$.

###### Proof.

Let $u_{0}\in B_{L^{\infty}}(R)$ and $k\geq k_{\epsilon}$.
We first prove (i).
By Hölder’s inequality and the second condition in (23),
for any $u\in B_{L^{\infty}(0,T;L^{\infty})}(M^{\prime})$, we estimate

$$ $\displaystyle\|\Phi_{N_{k}}[u]\|_{L^{\infty}_{t}(0,T;L^{\infty}_{x})}$ $\displaystyle\leq\left\|\|G_{N_{k}}(t,x,y)\|_{L^{1}_{y}}\right\|_{L_{t}^{ \infty}(0,T;L_{x}^{\infty})}\|u_{0}\|_{L^{\infty}}$ $\displaystyle\qquad+\left\|\|G_{N_{k}}(t-\tau,x,y)\|_{L^{1}_{\tau}(0,t;L^{1}_{ y})}\right\|_{L_{t}^{\infty}(0,T;L_{x}^{\infty})}\|F(u)\|_{L^{\infty}(0,T;L^{ \infty})}$ $\displaystyle\leq 2C_{\mathcal{L}}R+2C_{\mathcal{L}}TC_{F}M^{\prime p}\leq M^{ \prime}.$ $$

Next, we prove (ii). By Hölder’s inequality and the third condition in (23), for
any $u\in B_{L^{\infty}(0,T;L^{\infty})}(M^{\prime})$, we have

$$ $\displaystyle\Bigl{|}\Phi[u](t,x)-\Phi_{N_{k}}[u](t,x)\Bigr{|}$ $\displaystyle\leq\|G(t,x,y)-G_{N_{k}}(t,x,y)\|_{L^{s^{\prime}}_{y}}\|u_{0}\|_{ L^{s}}$ $\displaystyle~{}~{}+\|G(t-\tau,x,y)-G_{N_{k}}(t-\tau,x,y)\|_{L^{r^{\prime}}_{ \tau}(0,t;L^{s^{\prime}}_{y})}\|F(u)\|_{L^{r}(0,T;L^{s})},$ $$

which implies that

$$ $\displaystyle\left\|\Phi[u]-\Phi_{N_{k}}[u]\right\|_{L^{r}(0,T;L^{s})}$ $\displaystyle\leq\left\|\|G(t,x,y)-G_{N}(t,x,y)\|_{L^{s^{\prime}}_{y}}\right\| _{L^{r}_{t}(0,T;L^{s}_{x})}\|u_{0}\|_{L^{s}}$ $\displaystyle\quad+\left\|\|G(t-\tau,x,y)-G_{N}(t-\tau,x,y)\|_{L^{r^{\prime}}_ {\tau}(0,t;L^{s^{\prime}}_{y})}\right\|_{L^{r}_{t}(0,T;L^{s}_{x})}\|F(u)\|_{L^ {r}(0,T;L^{s})}$ $\displaystyle\leq C^{\prime}_{G}(N_{k})\|u_{0}\|_{L^{s}}+C_{G}(N_{k})\|F(u)\|_ {L^{r}(0,T;L^{s})}$ $\displaystyle\leq\epsilon|D|^{\frac{1}{s}}\|u_{0}\|_{L^{\infty}}+\epsilon\|F(u )\|_{L^{r}(0,T;L^{s})}$ $\displaystyle\leq\left(|D|^{\frac{1}{s}}R+C_{F}M^{p}\right)\epsilon$ $\displaystyle=C_{1}\epsilon,$ $$

where we note that $B_{L^{\infty}(0,T;L^{\infty})}(M^{\prime})\subset B_{L^{r}(0,T;L^{s})}(M)$.
Therefore, (ii) is proved. Thus, the proof of Lemma 2 is finished.
∎

Step 2.
We define the map $\Phi_{N_{k},net}$ by

$$ $\Phi_{N_{k},net}[u](t,x):=\int_{D}G_{N_{k}}(t,x,y)u_{0}(y)dy+\int_{0}^{t}\int_ {D}G_{N_{k}}(t-\tau,x,y)F_{net}(u(\tau,y))d\tau dy,$ (24) $$

where $F_{net}:\mathbb{R}\to\mathbb{R}$ is a ReLU neural network.
Here, as $F|_{(-M^{\prime},M^{\prime})}\in W^{1,\infty}(-M^{\prime},M^{\prime})$ from Assumption 2,
we can apply the approximation result by
[, Theorem 4.1] to see that
there exists
a ReLU neural network $F_{net}:\mathbb{R}\to\mathbb{R}$
such that

$$ $\|F-F_{net}\|_{L^{\infty}(-M,M)}\leq\epsilon,$ (25) $$

where the depth $L=L(F_{net})$ and the number $H=H(F_{net})$ of neurons are estimated as

$$ $\begin{cases}L(F_{net})\leq C\log_{2}\left(\epsilon^{-1}\right),\\ H(F_{net})\leq C\epsilon^{-1}\log_{2}\left(\epsilon^{-1}\right),\end{cases}$ (26) $$

where $C=C(M^{\prime},F)>0$ is a constant depending on $M^{\prime},F$.

###### Lemma 3 .

For any $u_{0}\in B_{L^{\infty}}(R)$ and $k\geq k_{\epsilon}$, the following statements hold:

- (i)
$\Phi_{N_{k},net}$ is a map from $B_{L^{\infty}(0,T;L^{\infty})}(M^{\prime})$ to itself.
- (ii)
There exists a constant $C_{2}>0$ such that
$\|\Phi_{N_{k}}[u]-\Phi_{N_{k},net}[u]\|_{L^{r}(0,T;L^{s})}\leq C_{2}\epsilon.$
for $u\in B_{L^{\infty}(0,T;L^{\infty})}(M^{\prime})$. Here, $C_{2}=2C_{\mathcal{L}}|D|^{\frac{1}{s}}T^{1+\frac{1}{r}}$.

###### Proof.

The proof follows the same lines with Lemma 2.
Let $u_{0}\in B_{L^{\infty}}(R)$ and $k\geq k_{\epsilon}$.
By Hölder’s inequality and Assumption 2,
for any $u\in B_{L^{\infty}(0,T;L^{\infty})}(M^{\prime})$, we estimate

$$ $\displaystyle\|\Phi_{N_{k},net}[u]\|_{L^{\infty}(0,T;L^{\infty})}$ $\displaystyle\leq\left\|\|G_{N_{k}}(t,x,y)\|_{L^{1}_{y}}\right\|_{L^{\infty}_{ t}(0,t;L^{\infty}_{x})}\|u_{0}\|_{L^{\infty}}$ $\displaystyle~{}~{}+\left\|\|G_{N_{k}}(t-\tau,x,y)\|_{L^{1}_{\tau}(0,t;L^{1}_{ y})}\right\|_{L^{\infty}_{t}(0,t;L^{\infty}_{x})}\|F_{net}(u)\|_{L^{\infty}(0, T;L^{\infty})}$ $\displaystyle\leq 2C_{\mathcal{L}}R+2C_{\mathcal{L}}T(1+C_{F}M^{\prime p})\leq M ^{\prime},$ $$

where we used the inequality

$$ $\displaystyle\|F_{net}(u)\|_{L^{\infty}(0,T;L^{\infty})}$ $\displaystyle\leq\|F(u)-F_{net}(u)\|_{L^{\infty}(0,T;L^{\infty})}+\|F(u)\|_{L^ {\infty}(0,T;L^{\infty})}$ $\displaystyle\leq\epsilon+C_{F}M^{\prime p}\leq 1+C_{F}M^{\prime p}.$ $$

Hence, (i) is proved.

Next, we prove (ii).
By Hölder’s inequality,
for any $u\in B_{L^{\infty}(0,T;L^{\infty})}(M^{\prime})$,
we estimate

$$ $\begin{split}&\Bigl{|}\Phi_{N_{k}}[u](t,x)-\Phi_{N_{k},net}[u](t,x)\Bigr{|}\\ &\leq\|G_{N_{k}}(t-\tau,x,y)\|_{L^{1}_{\tau}(0,t;L^{1}_{y})}\|F(u)-F_{net}(u) \|_{L^{\infty}(0,T;L^{\infty})}\\ &\leq 2C_{\mathcal{L}}T\|F(u)-F_{net}(u)\|_{L^{\infty}(0,T;L^{\infty})}\leq 2C _{\mathcal{L}}T\epsilon.\end{split}$ $$

Hence, we obtain

$$ $\begin{split}&\|\Phi_{N_{k}}[u]-\Phi_{N_{k},net}[u]\|_{L^{r}(0,T;L^{s})}\leq 2 C_{\mathcal{L}}|D|^{\frac{1}{s}}T^{1+\frac{1}{r}}\epsilon=C_{2}\epsilon.\end{split}$ $$

Therefore, (ii) is proved. Thus, the proof of Lemma 3 is finished.
∎

Step 3.
We define $\Gamma:B_{L^{\infty}}(R)\to B_{L^{r}(0,T;L^{s})}(M)$ by

$$ $\Gamma(u_{0}):=\Phi_{N_{k},net}^{[J]}[0]\quad\text{for }u_{0}\in B_{L^{\infty} }(R).$ (27) $$

###### Lemma 4 .

Let $J=\lceil\frac{\log(1/\epsilon)}{\log(1/\delta)}\rceil\in\mathbb{N}$.
Then there exists a constant $C_{3}>0$
such that for any $u_{0}\in B_{L^{\infty}}(R)$

$$ $\|\Gamma^{+}(u_{0})-\Gamma(u_{0})\|_{L^{r}(0,T;L^{s})}\leq C_{3}\epsilon,$ $$

where $C_{3}=M+\frac{C_{1}+C_{2}}{1-\delta}$.

###### Proof.

By the triangle inequality, we have

$$ $\|\Gamma^{+}(u_{0})-\Gamma(u_{0})\|_{L^{r}(0,T;L^{s})}\leq\|\Gamma^{+}(u_{0})- \Phi^{[J]}[0]\|_{L^{r}(0,T;L^{s})}+\|\Phi^{[J]}[0]-\Gamma(u_{0})\|_{L^{r}(0,T; L^{s})}.$ $$

As to the first term,
since $\Phi:B_{L^{r}(0,T;L^{s})}(M)\to B_{L^{r}(0,T;L^{s})}(M)$ is $\delta$-contractive and $u=\Gamma^{+}(u_{0})$ is the fixed point of $\Phi$ from Proposition 2, we have

$$ $\displaystyle\|\Gamma^{+}(u_{0})-\Phi^{[J]}[0]\|_{L^{r}(0,T;L^{s})}$ $\displaystyle=\|\Phi^{[J]}[u]-\Phi^{[J]}[0]\|_{L^{r}(0,T;L^{s})}$ $\displaystyle\leq\delta^{J}\|u-0\|_{L^{r}(0,T;L^{s})}\leq M\delta^{J}\leq M\epsilon.$ (28) $$

As to the second term, we see that

$$ $\displaystyle\|\Phi^{[J]}[0]-\Gamma(u_{0})\|_{L^{r}(0,T;L^{s})}=\|\Phi^{[J]}[0 ]-\Phi_{N_{k},net}^{[J]}[0]\|_{L^{r}(0,T;L^{s})}$ $\displaystyle\leq\sum_{j=1}^{J}\left\|\left(\Phi^{[J-j+1]}\circ\Phi_{N_{k},net }^{[j-1]}\right)[0]-\left(\Phi^{[J-j]}\circ\Phi_{N_{k},net}^{[j]}\right)[0] \right\|_{L^{r}(0,T;L^{s})}$ (29) $\displaystyle\leq\sum_{j=1}^{J}\delta^{J-j}\left\|\left(\Phi\circ\Phi_{N_{k}, net}^{[j-1]}\right)[0]-\Phi_{N_{k},net}^{[j]}[0]\right\|_{L^{r}(0,T;L^{s})}$ (30) $\displaystyle=\sum_{j=1}^{J}\delta^{J-j}\|\Phi[u_{j-1,N}]-\Phi_{N_{k},net}[u_{ j-1,N}]\|_{L^{r}(0,T;L^{s})},$ (31) $$

where $u_{j,N}:=\Phi_{N_{k},net}^{[j]}[0]$ and $u_{j-1,N}\in B_{L^{\infty}(0,T;L^{\infty})}(M^{\prime})\subset B_{L^{r}(0,T;L^
{s})}(M)$ by the third condition in (23).
By Lemma 2 and Lemma 3, we see that

$$ $\displaystyle\|\Phi[u_{j,N}]-\Phi_{N_{k},net}[u_{j,N}]\|_{L^{r}(0,T;L^{s})},$ $\displaystyle\leq\|\Phi[u_{j-1,N}]-\Phi_{N_{k}}[u_{j-1,N}]\|_{L^{r}(0,T;L^{s}) }+\|\Phi_{N_{k}}[u_{j-1,N}]-\Phi_{N_{k},net}[u_{j-1,N}]\|_{L^{r}(0,T;L^{s})}$ $\displaystyle\leq(C_{1}+C_{2})\epsilon.$ $$

Hence,

$$ $\displaystyle\|\Phi^{[J]}[0]-\Gamma(u_{0})\|_{L^{r}(0,T;L^{s})}\leq\sum_{j=1}^ {J}\delta^{J-j}(C_{1}+C_{2})\epsilon\leq\sum_{j=0}^{\infty}\delta^{j}(C_{1}+C_ {2})\epsilon=\frac{C_{1}+C_{2}}{1-\delta}\epsilon.$ (32) $$

Therefore, by (28) and (32), we conclude that

$$ $\displaystyle\|\Gamma^{+}(u_{0})-\Gamma(u_{0})\|_{L^{r}(0,T;L^{s})}$ $\displaystyle\leq\|\Gamma^{+}(u_{0})-\Phi^{[J]}[0]\|_{L^{r}(0,T;L^{s})}+\|\Phi ^{[J]}[0]-\Gamma(u_{0})\|_{L^{r}(0,T;L^{s})}$ $\displaystyle\leq\left\{M+\frac{C_{1}+C_{2}}{1-\delta}\right\}\epsilon=C_{3}\epsilon.$ $$

Thus, the proof of Lemma 4 is finished.
∎

Final step.
Finally, it is sufficient to represent the approximate operator $\Gamma$ given in (27) as a neural operator in the form of Definition 1 and to provide its quantitative estimates.

###### Lemma 5 .

$$ $\Gamma\in\mathcal{NO}^{L,H,ReLU}_{N,\varphi,\psi},$ $$

where $L=L(\Gamma)$, $H=H(\Gamma)$ satisfies

$$ $L(\Gamma)\leq C(\log(\epsilon^{-1}))^{2}\text{~{}~{}and~{}~{}}H(\Gamma)\leq C \epsilon^{-1}(\log(\epsilon^{-1}))^{2},$ $$

where $C>0$ is a constant depending on $\nu,M,F,D,p,d,r,s,R,\mathcal{L}$, and $N=N(\Gamma)$ satisfies

$$ $C_{G}(N)\leq\epsilon\text{~{}~{}and~{}~{}}C^{\prime}_{G}(N)\leq\epsilon.$ $$

###### Proof.

Let $N\geq N_{k}$. Then we have

$$ $\displaystyle C_{G}(N)\leq C_{G}(N_{k})\leq\epsilon,$ $\displaystyle C^{\prime}_{G}(N)\leq C^{\prime}_{G}(N_{k})\leq\epsilon.$ $$

We see that

$$ $\displaystyle\Phi_{N_{k},net}[u](t,x)$ $\displaystyle=\int_{D}G_{N_{k}}(t,x,y)u_{0}(y)dy+\int_{0}^{t}\int_{D}G_{N_{k}} (t-\tau,x,y)F_{net}(u(\tau,y))d\tau dy$ $\displaystyle=\sum_{n,m\in\Lambda_{N_{k}}}c_{n,m}\langle\psi_{m}(0,\cdot),u_{0 }\rangle\varphi_{n}(t,x)+\sum_{n,m\in\Lambda_{N_{k}}}c_{n,m}\langle\psi_{m},F_ {net}(u)\rangle\varphi_{n}(t,x)$ $\displaystyle=\sum_{n,m\in\Lambda_{N}}\tilde{c}_{n,m}\langle\psi_{m}(0,\cdot), u_{0}\rangle\varphi_{n}(t,x)+\sum_{n,m\in\Lambda_{N}}\tilde{c}_{n,m}\langle \psi_{m},F_{net}(u)\rangle\varphi_{n}(t,x),$ $$

where

$$ $\tilde{c}_{n,m}:=\begin{dcases}c_{n,m}&\text{if }n,m\in\Lambda_{N_{k}},\\ 0&\text{otherwise}.\end{dcases}$ $$

We denote by

$$ $\tilde{u}_{1}(t,x):=\sum_{n,m\in\Lambda_{N}}\tilde{c}_{n,m}\langle\psi_{m}(0, \cdot),u_{0}\rangle\varphi_{n}(t,x).$ $$

Then we see that

$$ $\Gamma(u_{0})(t,x):=\Phi^{[J]}_{N_{k},net}[0](t,x)=v_{J}(t,x),$ $$

where $v_{0}:=0$ and

$$ $\displaystyle v_{j+1}(t,x):=\tilde{u}_{1}(t,x)+\sum_{n,m\in\Lambda_{N}}\tilde{ c}_{n,m}\langle\psi_{m},F_{net}(v_{j})\rangle\varphi_{n}(t,x),\quad j=0,\cdots ,J-1.$ $$

We define $\tilde{K}^{(0)}_{N}:L^{\infty}(D)\to L^{r}(0,T;L^{s}(D))^{2}$ by

$$ $(\tilde{K}^{(0)}_{N}u_{0})(t,x):=\sum_{n,m\in\Lambda_{N}}\tilde{C}_{n,m}^{(0)} \langle\psi_{m}(0,\cdot),u_{0}\rangle\varphi_{n}(t,x)=\begin{pmatrix}\tilde{u} _{1}(t,x)\\ \tilde{u}_{1}(t,x)\end{pmatrix},$ $$

where $\tilde{C}_{n,m}^{(0)}=\begin{pmatrix}\tilde{c}_{n,m}\\
\tilde{c}_{n,m}\end{pmatrix}\in\mathbb{R}^{2\times 1}$.
We also define $\tilde{b}^{(0)}_{N}\in L^{r}(0,T;L^{s}(D))^{2}$ by

$$ $\tilde{b}^{(0)}_{N}(t,x):=\begin{pmatrix}\displaystyle\sum_{n,m\in\Lambda_{N}} \tilde{c}_{n,m}^{(0)}\langle\psi_{m},F_{net}(0)\rangle\varphi_{n}(t,x)\\ 0\end{pmatrix}=\sum_{n\in\Lambda_{N}}\tilde{b}_{N,n}^{(0)}\varphi_{n}(t,x),$ $$

where

$$ $\tilde{b}_{N,n}^{(0)}=\begin{pmatrix}\displaystyle\tilde{c}_{n,m}^{(0)}\langle \psi_{m},F_{net}(0)\rangle\\ 0\end{pmatrix}\in\mathbb{R}^{2\times 1}.$ $$

Then, we see that

$$ $(\tilde{K}^{(0)}_{N}u_{0})(t,x)+\tilde{b}^{(0)}_{N}(t,x)=\begin{pmatrix}v_{1}( t,x)\\ \tilde{u}_{1}(t,x)\end{pmatrix},$ $$

which corresponds to the input layer.

Next, we define

$$ $\tilde{W}=\begin{pmatrix}0&1\\ 0&1\end{pmatrix}\in\mathbb{R}^{2\times 2}$ $$

and

$$ $(\tilde{K}_{N}u)(t,x):=\sum_{n,m\in\Lambda_{N}}\tilde{C}_{n,m}\langle\psi_{m}, u\rangle\varphi_{n}(t,x)=\begin{pmatrix}\displaystyle\sum_{n,m\in\Lambda_{N}} \tilde{c}_{n,m}\langle\psi_{m},u_{1}\rangle\varphi_{n}(t,x)\\ 0\end{pmatrix}$ $$

for $u=(u_{1},u_{2})\in L^{r}(0,T;L^{s}(D))^{2}$,
where $\tilde{C}_{n,m}=\begin{pmatrix}\tilde{c}_{n,m}&0\\
0&0\end{pmatrix}\in\mathbb{R}^{2\times 2}$.
We also define $\tilde{F}_{net}:\mathbb{R}^{2}\to\mathbb{R}^{2}$ by

$$ $\tilde{F}_{net}(u)=\begin{pmatrix}F_{net}(u_{1})\\ u_{2}\end{pmatrix},\quad u=(u_{1},u_{2})\in\mathbb{R}^{2},$ $$

which is a ReLU neural network with (26) (Note that ReLU neural networks represent the identity map).
Then we write

$$ $\displaystyle\left[(\tilde{W}+\tilde{K}_{N})\circ\tilde{F}_{net}\begin{pmatrix }v_{j}\\ \tilde{u}_{1}\end{pmatrix}\right](t,x)$ $\displaystyle=\tilde{W}\begin{pmatrix}F_{net}(v_{j})(t,x)\\ \tilde{u}_{1}(t,x)\end{pmatrix}+\tilde{K}_{N}\begin{pmatrix}F_{net}(v_{j})\\ \tilde{u}_{1}\end{pmatrix}(t,x)$ $\displaystyle=\begin{pmatrix}\tilde{u}_{1}(t,x)+\displaystyle\sum_{n,m\in \Lambda_{N}}\tilde{c}_{n,m}\langle\psi_{m},F_{net}(v_{j})\rangle\varphi_{n}(t, x)\\ \tilde{u}_{1}(t,x)\end{pmatrix}$ $\displaystyle=\begin{pmatrix}v_{j+1}(t,x)\\ \tilde{u}_{1}(t,x)\end{pmatrix}\,\text{ for }\,j=1,\ldots,J-1.$ $$

which corresponds to the hidden layer.

Finally, denoting by

$$ $\tilde{W}^{\prime}:=(1,0)\in\mathbb{R}^{1\times 2},$ $$

which corresponds to the last layer, we obtain

$$ $\displaystyle\Gamma(u_{0})=\tilde{W}^{\prime}\circ\left[\underbrace{(\tilde{W} +\tilde{K}_{N})\circ\tilde{F}_{net}\circ\cdots\circ(\tilde{W}+\tilde{K}_{N}) \circ\tilde{F}_{net}}_{J-1}\right]\circ(\tilde{K}^{(0)}_{N}+\tilde{b}^{(0)}_{N })(u_{0})$ $$

By the above construction, we can check that $\Gamma\in\mathcal{NO}^{L,H,ReLU}_{N,\varphi,\psi}$.
Moreover, the depth $L(\Gamma)$ and the number $H(\Gamma)$ of neurons of the neural operator $\Gamma$ can be estimated as

$$ $\displaystyle L(\Gamma)\lesssim J\cdot L(F_{net})\lesssim\log(\epsilon^{-1}) \cdot\log_{2}(\epsilon^{-1})\lesssim(\log(\epsilon^{-1}))^{2},$ $\displaystyle H(\Gamma)\lesssim J\cdot H(F_{net})\lesssim\log(\epsilon^{-1}) \cdot\epsilon^{-1}\log_{2}(\epsilon^{-1})\lesssim\epsilon^{-1}(\log(\epsilon^{ -1}))^{2}.$ $$

Thus, the proof of Theorem 1 is complete.
∎

###### Lemma 1 .

###### Proof.

###### Lemma 2 .

###### Proof.

###### Lemma 3 .

###### Proof.

###### Lemma 4 .

###### Proof.

###### Lemma 5 .

###### Proof.

## Appendix E Proof of Corollary 1

By the triangle inequality, we estimate

$$ $\|u_{\max}-\hat{u}_{\max}\|_{L^{r}(0,\kappa^{*}T;L^{s})}\leq\sum_{\kappa=1}^{ \kappa^{*}}\|{\Gamma^{+}}^{[\kappa]}(u_{0})-\Gamma^{[\kappa]}(u_{0})\|_{L^{r}( (\kappa-1)T,\kappa T;L^{s})}.$ $$

Therefore, as $\kappa^{*}$ is finite, it is enough to prove that

$$ $\begin{split}&\|{\Gamma^{+}}^{[\kappa]}(u_{0})-\Gamma^{[\kappa]}(u_{0})\|_{L^{ r}((\kappa-1)T,\kappa T;L^{s})}\\ &\leq C\epsilon+\|{\Gamma^{+}}^{[\kappa-1]}(u_{0})((\kappa-1)T)-\Gamma^{[ \kappa-1]}(u_{0})((\kappa-1)T)\|_{L^{\infty}}\end{split}$ (33) $$

for some $C>0$ and for $\kappa=1,\ldots,\kappa^{*}$.
For $\kappa=1$,
the estimate (33) is already obtained in Theorem 1.
For $\kappa=2,\ldots,\kappa^{*}$, we have

$$ $\begin{split}&\|{\Gamma^{+}}^{[\kappa]}(u_{0})-\Gamma^{[\kappa]}(u_{0})\|_{L^{ r}((\kappa-1)T,\kappa T;L^{s})}\\ &\leq\|\Gamma^{+}({\Gamma^{+}}^{[\kappa-1]}(u_{0}))-\Gamma^{+}(\Gamma^{[\kappa -1]}(u_{0}))\|_{L^{r}((\kappa-1)T,\kappa T;L^{s})}\\ &\qquad+\|\Gamma^{+}(\Gamma^{[\kappa-1]}(u_{0}))-\Gamma(\Gamma^{[\kappa-1]}(u_ {0}))\|_{L^{r}((\kappa-1)T,\kappa T;L^{s})}.\end{split}$ $$

By Corollary 2 (or Proposition 1), the first term can be estimated as

$$ $\begin{split}&\|\Gamma^{+}({\Gamma^{+}}^{[\kappa-1]}(u_{0}))-\Gamma^{+}(\Gamma ^{[\kappa-1]}(u_{0}))\|_{L^{r}((\kappa-1)T,\kappa T;L^{s})}\\ &\leq C\|{\Gamma^{+}}^{[\kappa-1]}(u_{0})((\kappa-1)T)-\Gamma^{[\kappa-1]}(u_{ 0})((\kappa-1)T)\|_{L^{\infty}}.\end{split}$ $$

By Theorem 1, the second term in RHS can be estimated by $C\epsilon$ since $\Gamma^{[\kappa-1]}(u_{0})\in B_{L^{\infty}}(R)$.
Thus, the estimates (33) are proved for all $\kappa=1,\ldots,\kappa^{*}$.
∎

## Appendix F Preserving the positivity

We define the ReQU activation function by

$$ $ReQU(t):=\max\{0,t\}^{2},\quad t\in\mathbb{R}.$ $$

###### Corollary 5 .

Suppose that $\mathcal{L}$ and $F$ satisfy Assumptions 1, and 2, respectively, and the parameters $r,s$ satisfies the condition (3).
Assume that $F\in C^{1}(\mathbb{R};\mathbb{R})$ is a polynomial and satisfies that

$$ $F(t)\lesseqgtr 0\quad\text{ if $t\lesseqgtr 0$}.$ (34) $$

Moreover, assume that there exists $N_{0}\in\mathbb{N}$ such that the truncated expansion $G_{N}$ defined in Assumption 3 of the Green function $G$ satisfies

$$ $G_{N}(t,x,y)\geq 0,\quad 0<t<T,\ x,y\in D,$ $$

for any $N\geq N_{0}$.
Then, there exists $\Gamma\in\mathcal{NO}^{L,H,ReQU}_{N,\varphi,\psi}$ that the statement in Theorem 1 holds such that $\Gamma$ preserve the positivity
, i.e.,

$$ $\Gamma(u_{0})\lesseqgtr 0\quad\text{ if $u_{0}\in L^{\infty}(D)$ and $u_{0} \lesseqgtr 0$}.$ $$

###### Proof.

We can show that there exist $\Gamma\in\mathcal{NO}^{L,H,ReQU}_{N,\varphi,\psi}$ satisfying the statement of Lemma 4 where $G_{N_{k}}\geq 0$ and $F_{net}:\mathbb{R}\to\mathbb{R}$ is ReQU neural network and exactly represents the polynomial $F$ (see e.g., [, Theorem 2.2]).
Note that $\Gamma(u_{0})=\Phi_{N_{k},net}^{[J]}[0]$ where $\Phi_{N_{k},net}$ is a map defined in (24) and depends on $u_{0}$.
We will show that by induction

$$ $\Phi_{N_{k},net}^{[J]}[0]\lesseqgtr 0\ \text{ if }u_{0}\lesseqgtr 0.$ $$

Let $u_{0}\lesseqgtr 0$.
First, we see that $\Phi_{N_{k},net}[0]=\Phi_{N_{k},net}^{[1]}[0]\lesseqgtr 0$. Assume that $\Phi_{N_{k},net}^{[j]}[0]\lesseqgtr 0$ ($1\leq j\leq J-1$).
Since $F=F_{net}$ satisfies (34), we see that

$$ $\displaystyle\Phi_{N_{k},net}^{[j+1]}[0]$ $\displaystyle=\Phi_{N_{k},net}\left[\Phi_{N_{k},net}^{[j]}[0]\right]$ $\displaystyle=\underbrace{\int_{D}G_{N_{k}}(t,x,y)u_{0}(y)dy}_{\lesseqgtr 0}+ \underbrace{\int_{0}^{t}\int_{D}G_{N_{k}}(t-\tau,x,y)F_{net}(\Phi_{N_{k},net}^ {[j]}[0](\tau,y))d\tau dy}_{\lesseqgtr 0}\lesseqgtr 0.$ $$

Thus, the required result is proved.
∎

###### Corollary 5 .

###### Proof.

## Peer Reviews

confidence: 5
contribution: 3
initial_rating: 5
presentation: 4
questions: 1. Can the authors provide numerical experiments to validate the theoretical results and compare the performance of their neural operator with traditional PDE solvers or other neural operator architectures?
2. The dual of $L^\infty$ is more complex than that of $L^1$. The authors should consider revising or removing this aspect throughout the paper.
3. Could the authors address how the proposed method handles periodic boundary conditions, which are essential for applying FNOs? Additionally, could they discuss the role of the trace operator and how to ensure sufficient regularity for the boundary conditions when defining the neural operator $\Gamma$?
4. Can the authors elaborate on how the proposed framework can be extended to handle more general PDEs, specifically the Navier-Stokes equations mentioned in the abstract with their specific projection operator and boundary conditions?
5. In PDE theory, extending local solutions to global solutions is a challenging and significant issue. Although the paper acknowledges the challenge of extending local solutions to global ones, it does not provide concrete solutions or insights. Could the authors elaborate on this challenge and clarify the practical relevance of Corollary 1, especially regarding its error estimate for long-time solutions?
6. How does the choice of basis functions $(\phi$ and $\psi)$ in the neural operator affect the approximation error rates? Could the authors provide specific examples to illustrate these rates for different basis functions?
7. Regularity of the solution space is crucial in numerical approximation. However, this paper does not incorporate this aspect into an error analysis. Could the authors characterize the impact of function regularity on the error estimate?
soundness: 3
strengths: 1. The paper presents a quantitative approximation theorem for approximating solution operators of nonlinear parabolic PDEs using neural operators.
2. The authors show that the depth and number of neurons in their neural operators do not grow exponentially with the desired accuracy, potentially mitigating the "curse of parametric complexity."
3. The proof is constructive, leveraging the connection between Picard's iteration and the forward propagation through the layers of the neural operator.
4. The proposed framework may be potentially generalizable to other nonlinear PDEs solvable by Picard's iteration, including the Navier-Stokes, nonlinear Schrodinger, and nonlinear wave equations.
summary: This research focuses on a new method for approximating the solution operators of nonlinear parabolic PDEs using neural operators. The authors aim to bridge the gap between theoretical understanding of neural operators and their practical application as PDE solvers by developing a quantitative approximation theorem. This theorem demonstrates that neural operators, specifically designed with a structure inspired by Picard's iteration, can approximate solutions of parabolic PDEs efficiently without experiencing the exponential growth in model complexity often associated with general operator learning. The authors provide examples of how their method can be applied to Fourier neural operators (FNOs) and wavelet neural operators (WNOs), highlighting the potential for this approach to solve a wide range of nonlinear PDEs.
weaknesses: 1. The theoretical results are not supported by numerical experiments, making it difficult to assess the practical performance and efficiency of the proposed neural operator architecture.
2. The paper focuses on a specific class of parabolic PDEs amenable to analysis via Banach's fixed point theorem. The applicability to more general PDEs, particularly the Navier-Stokes equation mentioned in the abstract, is not demonstrated.
3. The paper does not adequately address the challenge of handling periodic boundary conditions, a crucial aspect for applying FNOs.
4. The paper does not discuss the role of the trace operator and how to ensure sufficient regularity for the boundary conditions when defining the neural operator $\Gamma$.
5. The assumptions on the operator L and the nonlinearity F are quite strong, potentially limiting the scope of applicability.

---

confidence: 4
contribution: 4
initial_rating: 8
presentation: 4
questions: 1. I do not think it is easy to find a Green's function or design one, which I see as the main drawback of the paper. Perhaps it is not that hard for the linear parabolic PDE. (addressed in the Questions below). How can one go about finding these (for specific cases or in general) or is this mainly constructive? Or do you only need to characterize the regularity of the Green function and how can you do that, if that is the case?
2.  Could you clarify if the error being log ((1/eps)^{-1}) implies exponential dependency?
3 (a). How can I interpret K^{(l)}_N? W are weights, b are biases, but I am not sure what K (the rank) is (unless it is specific to the Neural Operator?)
(b) Can you clarify (or add in the paper to make it more self-contained (see point 1 in Weakness) what is the rank of a neural operator? I do not believe it is defined in your paper (I checked other sources to understand). 
4. It seems the first 2 sentences of your Discussion are contradictory: you mention it is a limitation that your method only applies to parabolic PDEs, but the next sentence says your approach is not restricted to parabolic PDEs. Perhaps you could clarify that it is not truly a limitation then, and is just a first step?
5. For the long-time solution, how reasonable is it that u(T) again falls into  Ball(L_\infty, R) ? Does that limit the expressiveness of the solution of the PDE?
soundness: 4
strengths: This is a highly rigorous and technical paper. The details are presented clearly and I am able to understand the authors' motivation and thought process in developing the theory. I think this is a great theoretical contribution that has solid foundations in PDE theory to better understand Neural Operators (of which there are varieties like Fourier, Wavelet, etc.) and unifies all of them in this fairly decent class of PDEs. I hope to see extensions of this work (that the authors mentioned there would be) to more general classes of PDEs but this is a significant start. Using neural operators as solvers for PDEs has a wide array of applications, particularly in scientific computing and other "AI4Science" domains. What is interesting is the method and analysis they use to bypass, by appropriately selecting the basis functions within the neural operator, the exponential growth in model complexity. They also do not seem to have too many unreasonable assumptions (this I'm not highly confident about).  Error rates with respect to truncation N make sense. There is a lot of scope for future explorations with different architectures and PDEs.
summary: The authors derive the approximation rate of solution operators for the class of nonlinear parabolic partial differential equations (PDEs), contributing to the quantitative approximation theorem for solution operators of nonlinear PDEs. They show that neural operators can efficiently approximate these solution operators without the exponential growth in model complexity, thus strengthening the theoretical foundation of neural operators. An innovative link between Picard iteration and neural operators is found and the authors also leverage solid PDE theory (using Duhamel's principle). This work provides a very good theoretical foundation for neural operators (and generalizes previous works on Neural operators such as Fourier Neural Operators, etc.).
weaknesses: 1. The paper could provide a bit more intuition and background on Neural Operators to be a bit more self-contained. This could be presented in the Supplementary.
2. There are no numerical experiments (but for this paper I do not think this is a big weakness given that I believe the main contribution of theory is solid). Nevertheless, it would be nice to see in future (if it is even possible to implement all these different basis expansions, etc.)
3. I do not think it is easy to find a Green's function or design one, which I see as the main drawback of the paper (if there is to be a practical extension for Deep Learning). Perhaps it is not that hard to find one for the linear parabolic PDE. (addressed in the Questions below).
4. Minor: typo on line 182 "elliptic operato,r " 
I cannot think of any more now, but it depends on my questions below as well as I may not have understood everything.

---

confidence: 4
contribution: 3
initial_rating: 5
presentation: 3
questions: ### Additional points requiring clarification
1.	The paper mentions that Theorem 1 avoids the curse of dimensionality. However, this is not entirely correct: it seems that in the statement of the theorem, the involved constant $C$ may still grow exponentially in $d$?
2.	Theorem 1 requires $N$ to be chosen in such a way that $C_G(N) \leq \varepsilon$. Is it possible to quantitatively relate $N$ and $\varepsilon$, e.g., in the FNO case?
3.	P.3: You state that under suitable conditions on $\mathcal{L}$ and $F$, the problems (P) and (P’) are equivalent, if the function $u$ is sufficiently smooth. Could you give some precise references here?
4.	Appendix A does not fully clarify why the considered operators satisfy Assumption 1. For example, A.1 states a bound on the Green function, but does not make explicit why this bound implies Assumption 1. For readers less familiar with PDE theory it would be useful to make explicit (at least for one of the sections A.1, A.2, …) why this bound implies that Assumption 1 is satisfied. 
5.	Similar results to Proposition 1 are well-known in the literature. Can you relate Proposition 1 to the literature more precisely, outlining differences to existing results in the literature? 
6.	Fourier neural operators: why do you only consider the case $d=1$ here? 
7.	Could you clarify if WNOs with Haar wavelets, as defined in Appendix C.2, have been used in previous papers?
8.	What is the role of $\tilde{D}$ in Appendix C.1? Can’t we just set $\tilde{D}=D$?


### Additional comments:

-P.4, l182: it should be operator, not operato,r
soundness: 3
strengths: Overall, the paper is well-written and structured. The results hold in quite general settings: Assumption 1 on the differential operator $\mathcal{L}$ of the PDE is rather general and includes, e.g., Laplace operators with different types of boundary conditions. Assumption 2 on non-linearity seems more restrictive and mainly limited to power-type non-linearities, but the obtained results are useful even in the case $F=0$, which is also covered by the setting. Assumption 3 is crucial to make the approach work, based on an expansion of the underlying Green's function. In the supplementary material the paper specializes the setting to certain types of Fourier neural operators and Haar wavelet neural operators, demonstrating the usefulness of the result.
summary: The paper defines a neural operator based on Picard iteration and provides a quantitative approximation theorem for learning the solution operator of a general class of non-linear parabolic PDEs. The paper starts with a local well-posedness result for the considered PDEs, then introduces the considered neural operators, states the quantitative approximation theorem, provides a sketch of the proof and discusses extensions to longer time-horizons. The supplementary material contains proofs and shows how the result can be applied to Fourier Neural Operators (FNOs) and Haar Wavelet Neural Operators (WNOs).
weaknesses: The setting is formulated very abstractly, making reading potentially difficult for readers less familiar with general PDE theory. The paper does not discuss how the constructed neural operators could be implemented. In relation to existing neural operators (FNOs, WNOs) it remains unclear in what generality they are covered. Relation to existing literature is not always clear. I make these points more specific in the questions posed below.

Generally, I think this could be a valuable contribution, enhancing our understanding of mathematical foundations of neural operators. However, currently there seems to be a mistake in the proof of Theorem 1, affecting the construction of the neural operators and potentially making implementation infeasible. Consequently, at this point I need to recommend rejection, but I would be willing to increase my rating if this issue can be addressed convincingly. 

### Problem in the proof of Lemma 5

The main issue I currently see is that there appears to be a problem in the proof / formulation of the neural operators. The issue concerns the use of the time horizon $T$ instead of time variable $t$. In the definition of $\Phi_N$ in lines 264-269, in the first line the time integral is from $0$ to $t$ (the variable of the function to be approximated). In the next line, the integral with respect to $\tau$ is replaced by the inner product $\langle \psi_m, F(u) \rangle$. However, this means you are taking the integral up to $T$, not just up to $t$. 
The same problem can be found in the proof of Lemma 5 in the Appendix (lines 1245-1248), which does not seem to be correct in its current form. It seems that, to address this, the neural operator architecture would need to be modified, by making also the inner products $\langle \psi_m, u \rangle$ in the definition of the operator $K^{(\ell)}_N$ depend on the time variable $t$.
However, this would make a crucial difference from the perspective of computational effort required to evaluate the neural operator. In the current notation, to evaluate  $K^{(\ell)}_N$  only one inner product needs to be evaluated (the same for all times $t$). It seems that to fix the issue mentioned above, for each time $t$ a different inner product would need to be evaluated, which may be computationally infeasible (as at any given hidden layer, you are computing inner products with output functions from previous hidden layers). Here a thorough discussion is required.
Even without this issue, a more thorough discussion on implementation would be needed: The approach crucially relies on evaluating inner products, by definition of  $K^{(\ell)}_N[u]$ . To compute such inner products (between the hidden layer functions and the basis functions $\psi_m$) in practice, additional discretization is required at each layer; for each function evaluation.

---

confidence: 3
contribution: 3
initial_rating: 5
presentation: 3
questions: My primary concern is with the ReLU network approximation, specifically in the correspondence between the neural operator structure defined on page 6 and Lemma 2. In Eq. (24), you define the mapping, but it remains unclear how this mapping aligns with the neural operator structure described earlier. Please explain how the neural operator structure in page 6 result in eq (24). A clearer explanation or justification is needed to ensure consistency between these sections.

Additionally, I noticed that the notation $\Phi$ represents one step (either exact or approximated) of the iteration. However, on page 6, within each iteration, the activation is only used once, which is a shallow neural network with a single hidden layer. The reference to Guhring et al. (2020, Theorem 4.1) pertains to deep neural networks, which might not directly apply to this context. I suggest finding a more appropriate reference for the approximation result, specifically one that addresses shallow networks.
Another minor question is that, I think your remark “the curse of dimensionality that occurs in conventional neural network approximations does not appear here” in page 8 does not make sense. As addressed by yourself, the curse of dimensionality is hidden in Assumption 3, the representation of the Green’s function of the linear equation.

Another minor question is that, I think your remark “the curse of dimensionality that occurs in conventional neural network approximations does not appear here” in page 8 does not make sense. As addressed by yourself, the curse of dimensionality is hidden in Assumption 3, the representation of the Green’s function of the linear equation.
soundness: 3
strengths: Same as summary.
summary: This paper studies the approximation property of the neural operator for representing the solution operator for nonlinear parabolic equations. The key idea is that, the solution of the equation, written as the Duhamel’s integral, can be obtained from Picard iteration, with exponential convergence. The authors use a deep neural operator to approximate the Picard iteration process, where each layer approximates one step of Picard iteration. Overall, the paper is well written with clear explanation for all the ideas.
weaknesses: Need for clarification on ReLU approximation. Detail in Questions.

Meta-review: This work studies the approximation properties of neural operators quantitatively for a broad class of nonlinear parabolic equations. The main results demonstrate that a class of neural operators, inspired by Picard's iteration, can approximate solution operators of parabolic PDEs efficiently without the exponential growth in model complexity. The authors provide examples of how their framework can be applied to Fourier neural operators and wavelet neural operators, making a concrete connection between theoretical findings and practice in operator learning.

Overall, the paper is well-organized and clearly written. The assumptions and their generality/limitations are carefully explained, and the separation of results into local time and long-time settings enhances the logical flow of the arguments. Overall, this work makes a strong theoretical contribution, rooted in PDE theory, while offering a fresh perspective on operator learning. As such, the AC believes this paper warrants acceptance as a rigorous theoretical study, even in the absence of numerical experiments.

Decision: Accept (Poster)