# Padding-free Convolution based on Preservation of
Differential Characteristics of Kernels

Kuangdai Leng  Scientific Computing Department
Science and Technology Facilities Council
Didcot, UK
<EMAIL_ADDRESS>Jeyan Thiyagalingam  Scientific Computing Department
Science and Technology Facilities Council
Didcot, UK
<EMAIL_ADDRESS>
###### Abstract

Convolution is a fundamental operation in image processing and machine
learning. Aimed primarily at maintaining image size, padding is a key
ingredient of convolution, which, however, can introduce undesirable boundary
effects. We present a non-padding-based method for size-keeping convolution
based on the preservation of differential characteristics of kernels. The main
idea is to make convolution over an incomplete sliding window “collapse” to a
linear differential operator evaluated locally at its central pixel, which no
longer requires information from the neighbouring missing pixels. While the
underlying theory is rigorous, our final formula turns out to be simple: the
convolution over an incomplete window is achieved by convolving its nearest
complete window with a transformed kernel. This formula is computationally
lightweight, involving neither interpolation or extrapolation nor restrictions
on image and kernel sizes. Our method favours data with smooth boundaries,
such as high-resolution images and fields from physics. Our experiments
include: i) filtering analytical and non-analytical fields from computational
physics and, ii) training convolutional neural networks (CNNs) for the tasks
of image classification, semantic segmentation and super-resolution
reconstruction. In all these experiments, our method has exhibited visible
superiority over the compared ones.

###### Index Terms:

machine learning, computer vision, convolutional neural network, padding,
differential operator

##  I Introduction

Convolution is a basic operation in image processing. By convolving an image
with certain kernels, one can achieve various effects such as blurring,
sharpening, and edge detection. The establishment of modern convolutional
neural networks (CNNs) [1, 2, 3, 4] has unprecedentedly highlighted the
significance of convolution. CNN-based network architectures have been
thriving, such as variational autoencoders (VAEs) [5], generative adversarial
networks (GANs) [6] and more recently the diffusion models [7, 8], which spawn
numerous derivations for applications. Meanwhile, many techniques have been
developed to improve CNNs at a lower level, such as batch normalisation [9],
depth-wise separable convolution [10], and skip connections [2, 4], many of
which have become common practice in CNN tasks.

This paper concerns boundary handling, a key ingredient of convolution for
maintaining feature map size. Padding is the routine at present. Despite the
great success of CNNs with simple padding (e.g., zero padding), previous
studies have shown undesirable boundary effects caused by padding, such as
artefacts in features [11, 12, 13] and the spatial bias [14, 15]. A number of
techniques have been developed for alleviating padding-induced boundary
effects, such as explicit boundary handling (EBH) [11], partial convolution
[12], avoiding uneven application of padding [14], and quantifiable position-
information encoding [16]. Though they have been proven or supposed to improve
CNNs on different aspects, two major drawbacks persist. First, these
techniques only work with back propagation, hence unavailable for image
filtering with given kernels. Second, most of these techniques, along with
simple padding, are empirically motivated, lacking some rigorous connection
between the near-boundary and the interior parts whereby boundary handling can
be rendered more interpretable and controllable.

In this paper, we propose a padding-free method for size-keeping convolution,
available for both image filtering and CNN training. By introducing a
continuous presentation of image over a sliding window, we establish an
equivalence between window-wise convolution and pixel-wise differentiation.
The latter can be conducted locally at a boundary pixel so that padding is no
longer needed. Our final formula is elegantly simple and computationally
lightweight: the convolution over an incomplete window is achieved by
convolving its nearest complete window with a transformed kernel.

Boundary handling essentially addresses the issue of missing information, so
there cannot exist a method that always prevails for all kinds of data. The
reduction from convolution to differentiation makes our method more efficient
for images with smoother or more predictable boundaries, such as continuous
fields from mathematics and physics, and high-resolution images. Concerning
fields, a highly relevant application is CNN-based physics-informed learning
[17, 18, 19, 20, 21, 22], aimed at simulating or inverting partial
differential equations (PDEs) with physics-embedding loss functions. In a
physics-informed CNN, kernels are employed as a discrete representation of
differential operators, so preserving the differential characteristics of such
kernels at the boundary becomes imperative. First, any artefacts from padding
will act as secondary sources to the boundary-value problem, generating fake
energies to propagate across the entire domain; second, padding also means
imposing a Dirichlet boundary condition [22], which can be incompatible with
the given PDE system. Owing to its rich content, we have to discuss physics-
informed learning in another paper; here we propose our method for general-
purpose image filtering and machine learning.

The remainder of this paper is organised as follows. In the next section, we
cover the related work, followed by the description of our method in Section
III. We then describe our experiments in Section IV, covering both forward
filtering and CNN training. The paper is then concluded in Section V.

##  II Related Work

Most existing boundary handling techniques are motivated by CNNs, and thus
forward-incompatible. We found two methods available for forward image
filtering: padding by algebraic extrapolation [23] and a discrete Fourier
transformation-based method that involves padding by reflection and circular
deconvolution [24]. We will compare our method with the former, as we could
not find a reliable implementation of the latter.

For CNNs, studies have shown that simple padding may not only cause artefacts
in features [11, 12, 13] but also introduce a spatial bias that impairs the
translation invariance of CNNs [14, 15, 25, 26]. Here translation invariance
means that CNNs are expected to extract the relevant features regardless of
the absolute positions of entities in images. A metric of this bias has
recently been proposed in [16]. Existing remedies can be largely divided into
two categories: advanced padding and relative position encoding, the former
motivated more by artefact suppression and the latter by mitigating the
spatial bias. Advanced padding includes randomly-valued padding [27],
randomly-positioned padding [28], training an auxiliary CNN for padding [29],
symmetric padding with even-sized kernels [30], and some non-generic
algorithms for domain-specific data [31, 32, 13]. Relative position encoding
includes eliminating uneven application of padding by constraining the image
and kernel sizes [14], partial convolution [12, 33], and EBH [11]. Partial
convolution is a non-trainable method that first conducts convolution with
zero padding and then divides the result by the number fraction of existent
pixels in the sliding window. Statistically, it is similar to the randomly-
valued padding [27] where the padded values are sampled from a probability
distribution determined by a selected boundary vicinity. EBH is the most
expensive method, which introduces $(K^{2}-1)$ duplicates of kernels (where
$K$ is the kernel size) to be trained exclusively on the near-boundary pixels
grouped by their positions relative to the image boundary. In theory, these
$(K^{2}-1)$ duplicated kernels should maximally reduce boundary effects, but
such an extra cost is prohibitively high, especially for a large kernel size
(such as $K=7$); besides, these duplicated kernels can only see a small
fraction of data near the boundary, so they can converge much slower than the
main kernel trained for the bulk interior. Empirically, we do not observe an
outstanding advantage of EBH from our CNN experiments.

In Section IV, we will compare our method to the simple padding schemes
(zeros, reflect, replicate and circular), along with padding by extrapolation
[23], padding by distribution [27], partial convolution [12], and EBH [11].

##  III Method

We describe our method in this section. Einstein summation convention is
adopted for both superscript and subscript indices (in lower case letters)
unless they are parenthesised.

###  III-A Forward convolution

Let $\bm{\upomega}=\\{\omega_{ij}$}, $i,j\in\\{0,1,\cdots,K-1\\}$ denote the
input kernel, which has size $K\times K$, with $K$ being an odd number. Here
we assume a square-shaped kernel only to simplify the notations. Consider a
pixel “$\mathtt{a}$” centring a complete sliding window “$\mathtt{A}$”, as
illustrated in Fig. 1a, meaning that pixel $\mathtt{a}$ is “valid” for
convolution. Let $\mathbf{u}^{\mathtt{A}}=\\{u_{ij}^{\mathtt{A}}\\}$,
$i,j\in\\{0,1,\cdots,K-1\\}$ denote the input image given at the $K\times K$
pixels in $\mathtt{A}$ (after stride and dilation if required). The
convolution over $\mathtt{A}$ can then be written as

$z^{\mathtt{a}}:=\omega_{ij}u_{ij}^{\mathtt{A}}=\bm{\upomega}:\mathbf{u}^{\mathtt{A}}.$
(1)

We aim for a non-padding-based method to accomplish such convolution at the
“invalid” pixels that centre incomplete sliding windows, such as
“$\mathtt{b}$” centring “$\mathtt{B}$” in Fig. 1a.

A continuous, sub-pixel image can be formed in $\mathtt{A}$ using the Lagrange
interpolating polynomial, as denoted by $\tilde{u}^{\mathtt{A}}(h,w)$ with $h$
and $w$ being the spatial coordinates:

$\tilde{u}^{\mathtt{A}}(h,w):=l_{i}(h)l_{j}(w)u_{ij}^{\mathtt{A}},\quad
h,w\in[0,K-1],$ (2)

where $l_{i}(x)$ is the Lagrange basis simplified for a uniform grid with a
unit interval,

$l_{i}(x)=\prod_{\begin{subarray}{c}0\leq k\leq K-1\\\ k\neq
i\end{subarray}}\frac{x-k}{i-k},\quad i\in\\{0,1,\cdots,K-1\\}.$ (3)

The interpolation in eq. (2) leads to a 2D polynomial of degree $(K-1)$ that
preserves the pixel values, i.e.,
$\tilde{u}^{\mathtt{A}}(i,j)=u^{\mathtt{A}}_{ij}$. Spatial derivatives of
$\tilde{u}^{\mathtt{A}}(h,w)$ can then be conducted up to order $(K-1)$ in
each direction, formulated as

$\left.\mathcal{D}^{mn}\tilde{u}^{\mathtt{A}}\right|_{(h,w)}:=\left.\frac{\partial^{(m+n)}\tilde{u}^{\mathtt{A}}}{\partial
h^{(m)}\partial
w^{(n)}}\right|_{(h,w)}=l_{i}^{m}(h)l_{j}^{n}(w)u_{ij}^{\mathtt{A}},$ (4)

for $m,n\in\\{0,1,\cdots,K-1\\}$, where $l_{i}^{m}(x)$ is the $m$-th
derivative of $l_{i}(x)$, available in exact form given $K$ (as they are all
rational numbers), and $\mathcal{D}^{mn}$ denotes the partial differential
operator of order $(m,n)$, e.g.,
$\mathcal{D}^{12}=\frac{\partial^{3}}{\partial h\partial w^{2}}$. Evaluated at
a pixel with location $(p,q)$, for $p,q\in\\{0,1,\cdots,K-1\\}$, eq. (4)
yields

$\left.\mathcal{D}^{mn}\tilde{u}^{\mathtt{A}}\right|_{(p,q)}=\delta_{pq,ij}^{mn}u_{ij}^{\mathtt{A}}=\bm{\updelta}_{pq}^{mn}:\mathbf{u}^{\mathtt{A}},$
(5)

based on the definition that

$\delta_{pq,ij}^{mn}:=l_{i}^{m}(p)l_{j}^{n}(q).$ (6)

Equation (5) states that,
$\left.\mathcal{D}^{mn}\tilde{u}^{\mathtt{A}}\right|_{(p,q)}$, the $(m,n)$-th
order derivative of our continuous image evaluated at pixel $(p,q)$, can be
obtained by convolving the image $\mathbf{u}^{\mathtt{A}}$ with the above-
defined kernel $\bm{\updelta}^{mn}_{pq}$. Therefore, we call $\bm{\updelta}$
defined by eq. (6) the _differential kernels_ , which depends only on the
kernel size $K$. An example for $K=3$ is provided in Fig. 1b.

$u_{11}^{\mathtt{A}}$$u_{12}^{\mathtt{A}}$$u_{21}^{\mathtt{A}}$$u_{22}^{\mathtt{A}}$$u_{01}^{\mathtt{A}}$$u_{02}^{\mathtt{A}}$$u_{20}^{\mathtt{A}}$$u_{10}^{\mathtt{A}}$$u_{00}^{\mathtt{A}}$Pixel
$\mathtt{a}$Pixel $\mathtt{b}$Window $\mathtt{A}$Window $\mathtt{B}$(a) Image
setupPixels centring a complete sliding windowPixels centring an incomplete
sliding windowAbsent pixels outside image boundary(b) Differential kernels at
central pixel $\mathtt{a}$$\begin{aligned} &\begin{aligned}
&\bm{\updelta}^{00}_{11}=\begin{pmatrix}0&0&0\\\ 0&1&0\\\
0&0&0\end{pmatrix}&\Leftarrow&\quad\mathcal{D}^{00}=1\\\
&\bm{\updelta}^{01}_{11}=\dfrac{1}{2}\begin{pmatrix}0&0&0\\\ -1&0&1\\\
0&0&0\end{pmatrix}&\Leftarrow&\quad\mathcal{D}^{01}=\dfrac{\partial}{\partial
w}\\\ &\bm{\updelta}^{02}_{11}=\begin{pmatrix}0&0&0\\\ 1&-2&1\\\
0&0&0\end{pmatrix}&\Leftarrow&\quad\mathcal{D}^{02}=\dfrac{\partial^{2}}{\partial
w^{2}}\\\ &\bm{\updelta}^{11}_{11}=\dfrac{1}{4}\begin{pmatrix}1&0&-1\\\
0&0&0\\\
-1&0&1\end{pmatrix}&\Leftarrow&\quad\mathcal{D}^{11}=\dfrac{\partial}{\partial
h}\dfrac{\partial}{\partial w}\\\
&\bm{\updelta}^{12}_{11}=\dfrac{1}{2}\begin{pmatrix}-1&2&-1\\\ 0&0&0\\\
1&-2&1\end{pmatrix}&\Leftarrow&\quad\mathcal{D}^{12}=\dfrac{\partial}{\partial
h}\dfrac{\partial^{2}}{\partial w^{2}}\\\
&\bm{\updelta}^{22}_{11}=\begin{pmatrix}1&-2&1\\\ -2&4&-2\\\
1&-2&1\end{pmatrix}&\Leftarrow&\quad\mathcal{D}^{22}=\dfrac{\partial^{2}}{\partial
h^{2}}\dfrac{\partial^{2}}{\partial w^{2}}\\\ \end{aligned}\\\
&\bm{\updelta}^{10}_{11}=\left(\bm{\updelta}^{01}_{11}\right)^{\top},\quad\bm{\updelta}^{20}_{11}=\left(\bm{\updelta}^{02}_{11}\right)^{\top},\quad\bm{\updelta}^{21}_{11}=\left(\bm{\updelta}^{12}_{11}\right)^{\top}\end{aligned}$
(c) Assembled differential kernels: $\vec{\mathbf{D}}_{11}$ for central pixel
$\mathtt{a}$ and $\vec{\mathbf{D}}_{00}$ for target pixel $\mathtt{b}$
($R=S=0$)
$\vec{\mathbf{D}}_{11}=\begin{pmatrix}0&0&0&0&\frac{1}{4}&-\frac{1}{2}&0&-\frac{1}{2}&1\\\
0&0&0&-\frac{1}{2}&0&1&1&0&-2\\\
0&0&0&0&-\frac{1}{4}&-\frac{1}{2}&0&\frac{1}{2}&1\\\
0&-\frac{1}{2}&1&0&0&0&0&1&-2\\\ 1&0&-2&0&0&0&-2&0&-4\\\
0&\frac{1}{2}&1&0&0&0&0&-1&-2\\\
0&0&0&0&-\frac{1}{4}&\frac{1}{2}&0&-\frac{1}{2}&1\\\
0&0&0&\frac{1}{2}&0&-1&1&0&-2\\\
0&0&0&0&\frac{1}{4}&\frac{1}{2}&0&\frac{1}{2}&1\\\
\end{pmatrix},\quad\vec{\mathbf{D}}_{00}=\begin{pmatrix}1&-\frac{3}{2}&1&-\frac{3}{2}&\frac{9}{4}&-\frac{3}{2}&1&-\frac{3}{2}&1\\\
0&2&-2&0&-3&3&0&2&-2\\\
0&-\frac{1}{2}&1&0&\frac{3}{4}&-\frac{3}{2}&0&-\frac{1}{2}&1\\\
0&0&0&2&-3&2&-2&3&-2\\\ 0&0&0&0&4&-4&0&-4&4\\\ 0&0&0&0&-1&2&0&1&-2\\\
0&0&0&-\frac{1}{2}&\frac{3}{4}&-\frac{1}{2}&1&-\frac{3}{2}&1\\\
0&0&0&0&-1&1&0&2&-2\\\
0&0&0&0&\frac{1}{4}&-\frac{1}{2}&0&-\frac{1}{2}&1\end{pmatrix}$ (d) Input
kernel $\bm{\upomega}$ (a box-blur filter) and the transformed kernel
$\bm{\upvarpi}$$\bm{\upomega}=\begin{pmatrix}1&1&1\\\ 1&1&1\\\ 1&1&1\\\
\end{pmatrix}\quad\Longrightarrow\quad\bm{\upvarpi}=\begin{pmatrix}16&-8&4\\\
-8&4&-2\\\ 4&-2&1\\\ \end{pmatrix},\quad\text{by}\ \
\vec{\bm{\upvarpi}}=\vec{\mathbf{D}}_{00}\cdot\vec{\mathbf{D}}_{11}^{-1}\cdot\vec{\bm{\upomega}}$
Figure 1: A complete example of our method: convolving a $4\times 4$ image
with a $3\times 3$ kernel. We aim to conduct convolution at invalid pixel
$\mathtt{b}$, whose nearest complete window is $\mathtt{A}$. Based on eq.
(12), the outcome is given by
$z^{\mathtt{b}}=\bm{\upvarpi}:\mathbf{u}^{\mathtt{A}}$, where
$\vec{\bm{\upvarpi}}=\vec{\mathbf{D}}_{00}\cdot\vec{\mathbf{D}}_{11}^{-1}\cdot\vec{\bm{\upomega}}$,
as the location of $\mathtt{b}$ in $\mathtt{A}$ is $(R,S)=(0,0)$ and
$M=\frac{K-1}{2}=1$. The two shaded areas in (b) and (c) mark out that
$\bm{\updelta}^{02}_{11}$ constructs the third column of
$\vec{\mathbf{D}}_{11}$, according to eq. (10). The $\vec{\mathbf{D}}$ metrics
and their inverse are precomputed.

Our central idea is to _represent the wanted convolution or eq. ( 1) over a
generic window $\mathtt{B}$, complete or incomplete, as a unique linear
differential operator $\mathcal{L}$ applied on the continuous image
$\tilde{u}(h,w)$ and evaluated locally at the centre $\mathtt{b}$_. Formally,
we prescribe the following equivalence:

$\bm{\upomega}:\mathbf{u}^{\mathtt{B}}\equiv\left.\mathcal{L}\tilde{u}\right|_{\mathtt{b}},\quad\mathcal{L}:=\alpha^{mn}\mathcal{D}^{mn},$
(7)

where $\alpha^{mn}$, for $m,n\in\\{0,1,\cdots,K-1\\}$, are the real
coefficients of the linear differential operator $\mathcal{L}$. It must be
emphasised that eq. (7) does not specify how $\tilde{u}(h,w)$ is determined,
implying that the window for its interpolation does not need to be our target
window $\mathtt{B}$. This will eventually enable convolution over the
incomplete windows.

The coefficients $\alpha^{mn}$ are determined such that eq. (7) holds _at
every valid pixel_ given $\tilde{u}(h,w)$ interpolated by its centred window.
A generic example is our pixel $\mathtt{a}$ that centres window $\mathtt{A}$.
The local coordinates of $\mathtt{a}$ in $\mathtt{A}$ are $(M,M)$, with
$M=\frac{K-1}{2}$. With $\tilde{u}(h,w)$ interpolated by $\mathtt{A}$, eq. (7)
becomes

$\displaystyle\bm{\upomega}:\mathbf{u}^{\mathtt{A}}$
$\displaystyle\equiv\left.\mathcal{L}\tilde{u}\right|_{\mathtt{a}}\equiv\left.\mathcal{L}\tilde{u}^{\mathtt{A}}\right|_{\mathtt{a}}$
(8)
$\displaystyle=\left.\alpha^{mn}\mathcal{D}^{mn}\tilde{u}^{\mathtt{A}}\right|_{(M,M)}=\alpha^{mn}\bm{\updelta}_{MM}^{mn}:\mathbf{u}^{\mathtt{A}},$

the last part using our definition of $\bm{\updelta}^{mn}_{pq}$ in eq. (6). To
make eq. (8) hold regardless of data $\mathbf{u}^{\mathtt{A}}$, $\alpha^{mn}$
must be the solution of the following $K^{2}\times K^{2}$ linear system:

$\alpha^{mn}\bm{\updelta}_{MM}^{mn}=\bm{\upomega}.$ (9)

If we denote the vectorisation of a generic $K\times K$ matrix $\mathbf{a}$ by
$\vec{\mathbf{a}}$ such that $\vec{a}_{iK+j}=a_{ij}$, the above linear system
can be recast to the following standard form:

$\displaystyle\underbrace{\begin{pmatrix}\vec{\delta}_{MM,0}^{0}&\vec{\delta}_{MM,0}^{1}&\cdots&\vec{\delta}_{MM,0}^{K^{2}-1}&\\\
\vec{\delta}_{MM,1}^{0}&\vec{\delta}_{MM,1}^{1}&\cdots&\vec{\delta}_{MM,1}^{K^{2}-1}&\\\
\vdots&\vdots&\ddots&\vdots&\\\
\vec{\delta}_{MM,K^{2}-1}^{0}&\vec{\delta}_{MM,K^{2}-1}^{1}&\cdots&\vec{\delta}_{MM,K^{2}-1}^{K^{2}-1}&\end{pmatrix}}_{\begin{matrix}\vec{\mathbf{D}}_{MM}\end{matrix}}\underbrace{\begin{pmatrix}\vec{\alpha}^{0}\\\
\vec{\alpha}^{1}\\\ \vdots\\\
\vec{\alpha}^{K^{2}-1}\end{pmatrix}}_{\begin{matrix}\vec{\bm{\upalpha}}\end{matrix}}$
$\displaystyle\begin{aligned}
=\underbrace{\begin{pmatrix}\vec{\omega}_{0}&\vec{\omega}_{1}&\cdots&\vec{\omega}_{K^{2}-1}\end{pmatrix}^{\top}}_{\begin{matrix}\vec{\bm{\upomega}}^{\top}\end{matrix}}.\end{aligned}$
(10)

The assembled matrix $\vec{\mathbf{D}}_{MM}$ is always invertible because the
$K^{2}$ differential kernels, $\vec{\bm{\updelta}}^{k}_{MM}$ for
$k\in\\{0,1,\cdots,K^{2}-1\\}$ (i.e., each column in $\vec{\mathbf{D}}_{MM}$),
are linearly independent given that $l_{i}(x)$ is a complete polynomial of
degree $(K-1)$. The inverse of $\vec{\mathbf{D}}_{MM}$ can be exactly shown
for a given kernel size $K$, so computing
$\vec{\bm{\upalpha}}=\vec{\mathbf{D}}_{MM}^{-1}\cdot\vec{\bm{\upomega}}$ is
trivial.

Having uniquely determined the differential operator $\mathcal{L}$ as
$\alpha^{mn}\mathcal{D}^{mn}$, we can evaluate $\mathcal{L}$ at any invalid
pixel once a continuous image in its neighbourhood is provided. Here we choose
to determine such a continuous image by its nearest complete window. Let
$\mathtt{b}$ be an invalid pixel centring an incomplete window $\mathtt{B}$,
and $\mathtt{A}$ be the complete window nearest to $\mathtt{b}$, such as Fig.
1a. It is straightforward to show that $\mathtt{b}$ lies in $\mathtt{A}$ (but
not at its centre). Assume that the local coordinates of $\mathtt{b}$ in
$\mathtt{A}$ are $(R,S)$, with $R,S\in\\{0,1,\cdots,K-1\\}$ and
$(R,S)\neq(M,M)$. We use the R.H.S. of eq. (7) to compute the convolution over
$\mathtt{B}$, however, with the continuous image $\tilde{u}(h,w)$ interpolated
from $\mathtt{A}$:

$\displaystyle\bm{\upomega}:{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{\mathbf{u}^{\mathtt{B}}}}\equiv$
$\displaystyle\left.\mathcal{L}\tilde{u}\right|_{\mathtt{b}}\approx\left.\mathcal{L}\tilde{u}^{\mathtt{A}}\right|_{\mathtt{b}}$
(11) $\displaystyle=$
$\displaystyle\left.\alpha^{mn}\mathcal{D}^{mn}\tilde{u}^{\mathtt{A}}\right|_{(R,S)}=\alpha^{mn}\bm{\updelta}_{RS}^{mn}:\mathbf{u}^{\mathtt{A}}.$

We colour $\mathbf{u}^{\mathtt{B}}$ in red to indicate that it is partially
unavailable. Note that the “$\approx$” sign in the above equation indicates
the only approximation we have introduced: the window by which the continuous
image $\left.\tilde{u}\right|_{\mathtt{b}}$ is interpolated. More readably,
using
$\vec{\bm{\upalpha}}=\vec{\mathbf{D}}_{MM}^{-1}\cdot\vec{\bm{\upomega}}$, eq.
(11) can be simplified as

$\vec{\bm{\upomega}}\cdot{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\vec{\mathbf{u}}^{\mathtt{B}}}\approx\vec{\bm{\upvarpi}}\cdot\vec{\mathbf{u}}^{\mathtt{A}},\quad\vec{\bm{\upvarpi}}:=\vec{\mathbf{D}}_{RS}\cdot\vec{\mathbf{D}}_{MM}^{-1}\cdot\vec{\bm{\upomega}}.$
(12)

We refer to eq. (12) as the _differential kernel transformation_. Clearly, it
is compatible with the valid pixels, for which $R=S=M$ and thus
${\bm{\upvarpi}}={\bm{\upomega}}$. Refer to Fig. 1 for the example of $K=3$.
In summary, _the convolution with $\bm{\upomega}$ over an incomplete window
$\mathtt{B}$ is conducted by a “shifted” convolution with the transformed
kernel $\bm{\upvarpi}$ over its nearest complete window $\mathtt{A}$, under a
window shift of $(M-R,M-S)$_.

###  III-B Method properties

Our method has the following key properties:

####  III-B1 Theoretical soundness

Our method preserves the differential characteristics of kernels at the
invalid pixels by making the window-wise convolution collapse to a pixel-wise
differential operator, or eq. (7). Such a connection between the “invalid”
near-boundary part and the “valid” interior part is theoretically sound and
self-contained. In contrast, most previous methods are empirically motivated
by CNNs.

####  III-B2 Only using original pixel values

Though we have introduced a continuous image conceptually, our final formula
for convolution at the invalid pixels, eq. (12), operates on the original
pixels from the input image. Without introducing extra information outside the
image boundary (such as by padding or extrapolation) or between pixels (such
as by interpolation), our method avoids these sources of artefacts. This
property also makes our method compatible with stride and dilation for CNN
training.

####  III-B3 Low overhead

The transformation matrices in eq. (12),
$\vec{\mathbf{D}}_{RS}\cdot\vec{\mathbf{D}}_{MM}^{-1}$, are constants given
the kernel size $K$, which can be precomputed for frequently-used $K$’s, such
as 3, 5 and 7. Therefore, eq. (12) demands a low computational cost in both
forward and back propagation.

####  III-B4 Favouring data with smooth boundaries

Let $\tilde{u}^{*}(h,w)$ denote the true image function, which can be
discontinuous or even inexpressible. At the boundary pixels, our method
preserves the convolution-associated differential operator up to order
$(K-1)$. Therefore, it becomes exact when the local Taylor expansion of
$\tilde{u}^{*}(h,w)$ on the boundary is of order $(K-1)$ or lower, i.e., when
$\tilde{u}^{*}(h,w)$ is sufficiently smooth on the boundary. The error of our
method increases as $\tilde{u}^{*}(h,w)$ becomes more non-smooth or
unpredictable near the boundary. Datasets underpinned by physical or
mathematical processes, such as solutions of partial differential equations
and tomographic images, tend to benefit from our method. Furthermore, our
method tends to work better with higher-resolution images because
$\tilde{u}^{*}(h,w)$ becomes smoother as the sliding windows shrink with
respect to image contents.

##  IV Experiments

We evaluate our method with two types of experiments, the former on image
filtering with given kernels, as reported in Section IV-A, and the latter on
CNN-based computer vision tasks, as reported in Section IV-B. Einstein
summation convention is not used in this section.

###  IV-A Image filtering

We consider three synthetic datasets. The first two are analytical 2D
functions, respectively generated from the Chebyshev polynomials and the
spherical harmonics. Both are popular basis functions in computational physics
and mathematics. Therefore, the accuracy of our method as tested on these
basis functions can reasonably indicate its versatility for handling extensive
continuous fields. Our Chebyshev-based functions are given by

$C_{n}(h,w)=U_{n}(h)U_{n}(w)\sin\left(n(h+w)\right),$ (13)

where the $U_{n}$’s, for $n\in\\{0,1,2,\cdots\\}$, are the Chebyshev
polynomials of the second kind, and the rotation by $\sin\left(n(h+w)\right)$
makes the 2D patterns non-parallel to the axes. The spherical harmonic-based
functions are given by

$S_{n}(h,w)=Y_{2n}^{n}(h,w)\sin\left(n(h+w)\right),$ (14)

where $Y_{l}^{m}$ is the spherical harmonic function of degree $l$ and order
$m$, composed of trigonometric functions and the associated Legendre
polynomials. As the function order $n$ grows, both $C_{n}$ and $S_{n}$ becomes
more oscillating or non-smooth, as shown in Fig. 2a and 2b. Our third dataset
is non-analytical, including numerical solutions of the Navier–Stokes
equations for turbulence, as shown in 2c, borrowed from the physics-informed
neural operators [34].

We apply random filters to these datasets and compare the accuracy of the
following eight methods: padding respectively by zeros (Zero), reflection
(Refl), replicate (Repl), circular (Circ), extrapolation (Extr) [23], and
distribution (Rand) [27], along with partial convolution (Part) [12] and our
differentiation-based method (Diff). For Extr, we use linear, quadratic and
cubic respectively for $K=3,5$ and $7$. For Rand, the padded values are
sampled from the four normal distributions determined respectively for the
top, bottom, left and right edges with a thickness of $(K+1)/2$. Requiring
back-propagation, EBH [11] is not applicable here.

The superiority of our method against the others is visible from Fig. 2. The
second row shows the $L^{1}$ errors of each method for 100 random filters
applied to the three datasets. Our method proves to be remarkably more
accurate than the others (note that the y-axes are logarithmic). For Chebyshev
and spherical harmonics, the errors increase with the function order $n$, but
our method prevails across all the orders. The bottom row of Fig. 2 zooms into
the boundary artefacts caused by the Laplace filters with different kernel
sizes. It is shown that our method is visually artefact-free even at $n=100$;
Extr also works reasonably well, but its induced errors are still visible and
1$\sim$2 orders of magnitude larger than ours. The other padding schemes and
partial convolution will cause strong artefacts irrespective of kernel size
and function order, so they are unsuitable for the task of image filtering.

Figure 2: Filtering three datasets by convolution: Chebyshev, spherical
harmonics and solutions of Navier-Stokes equations [34]. The top row shows our
target functions. The second row shows the $L^{1}$ errors ($\varepsilon^{1}$)
for 100 random filters applied to the datasets. We vary the function order $n$
for Chebyshev and spherical harmonics, and the batch index for the Navier-
Stokes. In the bottom row, we apply the Laplace filters
($\Delta=\partial_{hh}+\partial_{ww}$) of sizes 3, 5 and 7 to $C_{100}$; (g)
shows the ground truth of $\Delta C_{100}$ based on analytical evaluation, and
(h)$\sim$(j) the results by convolution, zooming into the annotated box in
(g).

###  IV-B Learning with CNNs

In this section, we consider three common CNN tasks: image classification,
semantic segmentation and super-resolution reconstruction. For all these
experiments we use real-world datasets (instead of e.g. analytical fields) to
avoid favouritism towards our method by means of forward filtering. The
original U-Net architecture [35] is adopted, with its Conv2d layers varying
among nine different boundary handling methods. The first eight are those from
the previous experiment: Zero, Refl, Repl, Circ, Extr, Rand, Part and Diff
(ours), all non-trainable, and the last one is EBH [11], involving eight
duplicates of kernels (as $K=3$ in a U-Net). The relative wall-times for
training the nine U-Nets are reported in Table I, which shows that our method
runs as fast as circular padding (PyTorch built-in). For each problem, the
nine U-Nets are initialised with the same weights, and we loop over five
random seeds to obtain the reported metric scores. It must be emphasised that
_boundary handling is a low-level operation in CNNs, so we use a simple
network architecture and loss functions to isolate its influences_ instead of
pursuing the state of the art of the considered problems (datasets) with any
advanced yet irrelevant techniques.

####  IV-B1 Classification

We use the Caltech-101 dataset [36] for this experiment. The latent (bottom
layer) of the U-Net is connected to a two-layer fully-connected network to
predict the soft labels and then the classification (cross-entropy) loss. The
total loss is the sum of the classification loss and 1% of the reconstruction
loss (which accelerates convergence). The images are all reshaped to
$448\times 448$.

The accuracy of the models on the test set (20% of data) is reported in Table
I. It is shown that Part, EBH and Diff have achieved a much higher accuracy
($>50\%$) than the other six padding-based methods ($<30\%$), with our Diff
attaining the highest. The padding-based U-Nets have mostly failed to learn
(with many hyperparameters tested), as can be seen from their low accuracy and
training history (not shown here for brevity). This can be a good example of
padding-free boundary handling (Part, EBH and Diff) significantly enhancing
the learnability of a CNN.

TABLE I: Metric Scores for the CNN-based Experiments | Classification | Semantic Segmentation | Super-resolution Reconstruction | Computational cost
---|---|---|---|---
Method | Caltech-101 | Cityscapes | ETOPO-15′′ | Original U-Net
|  _Accuracy_ a | _IoU_ | _F1_ | _Accuracy_ | ${\varepsilon^{2}_{\text{INTER}}\times 10^{4}}$ | ${\varepsilon^{2}_{\text{FRAME}}\times 10^{4}}$ | ${\varepsilon^{2}_{\text{FRAME}}/\varepsilon^{2}_{\text{INTER}}}$ |  _# kernels_ b |  _Wall-time_ c
Zero | 30.7% | 84.8% | 91.8% | 97.9% | 7.0 | 10.2 | 144.5% | 1 | 1.0
Refl | 28.8% | 84.6% | 91.6% | 97.7% | 7.2 | 10.0 | 139.9% | 1 |  $\sim$1.1
Repl | 22.9% | 85.1% | 92.0% | 97.9% | 7.0 | _9.2_ | _132.0%_ | 1 |  $\sim$1.1
Circ | 29.1% | 85.0% | 91.9% | 97.9% | 7.2 | 11.6 | 157.9% | 1 |  $\sim$1.5
Extr [23] | 28.8% | 85.3% | 92.1% | 98.0% | 7.4 | 10.9 | 147.5% | 1 |  $\sim$5.8
Rand [27] | 27.9% | 85.1% | 91.9% | 97.9% | 7.7 | 15.5 | 200.2% | 1 |  $\sim$2.6
Part [12] | _53.7%_ | 85.5% | 92.3% | 98.1% | 7.1 | 10.3 | 145.0% | 1 |  $\sim$1.2
EBH [11] | 53.2% | _85.9%_ | _92.6%_ | _98.3%_ | _6.8_ | 9.3 | 135.9% | 9 |  $\sim$12.5
Diff (ours) |  55.8% | 86.1% | 92.8% | 98.4% | 6.6 | 8.0 | 120.7% | 1 |  $\sim$1.6
STD Diffd |  $\pm$0.92% |  $\pm$0.10% |  $\pm$0.06% |  $\pm$0.02% |  $\pm$0.05 |  $\pm$0.12 |  $\pm$1.8% |  |
a The best and second-best scores in each column are respectively printed in
boldface and italic-boldface.
b This column shows the total number of kernels to be trained for one target
kernel; only EBH introduces eight duplicates.
c This column shows the approximate wall-time (relative to Zero) required to
train a U-Net with input shape [64, 3, 224, 224], measured for forward
${}^{\mathrm{\ }}$ and back propagation on an Nvidia A100. The first four
methods are PyTorch built-ins while the rest are based on our implementation.
d For brevity, we only show the standard deviations of the metrics yielded by
Diff. The metrics are stable with respect to model initialisation.

####  IV-B2 Semantic Segmentation

In this experiment, we use the Cityscapes dataset [36] for end-to-end
supervised learning of semantic segmentation with a U-Net. The category
identities (8 classes) instead of the fine identities (34 classes) are used as
the labels because our simple architecture and loss function (cross entropy)
could not well handle a high degree of class imbalance in the latter. The
original images ($1024\times 2048$) are decimated by a factor of two due to
our device capacity.

The segmentation metric scores on the validation set are shown in Table I. It
can be seen that the scores from Zero, Refl, Repl, Circ and Rand are mostly
identical, implying that none of them have facilitated segmentation from the
perspective of boundary handling. Extr and Part have led to a small (yet
visible) improvement, and EBH has advanced further. Our method Diff has
yielded the best results. Note that all the methods have attained a high
baseline ($\text{accuracy}>97\%$) in an absolute sense, from which even a
small improvement is not easy to achieve.

####  IV-B3 Super-resolution Reconstruction

In this experiment, we train a U-Net to reconstruct a world topographic map
from a low to a high resolution. The data come from the _ETOPO 2022 15 Arc-
Second Global Relief Model_ [37], a large image containing $43200\times 86400$
pixels, each spanning a central angle of $15^{\prime\prime}$ (or 0.464 km on
Earth’s surface) in the latitudinal and longitudinal directions, as displayed
in Fig. 3a and 3b. We train the U-Net with non-overlapping small patches
sampled from this large image, each with size $192\times 192$ (geographically
$0.8^{\circ}\times 0.8^{\circ}$). We generate the low-resolution input by a
Gaussian filter ($\sigma=3$), from which we attempt to recover the high-
resolution output, as shown in Fig. 3c. Mean squared error (MSE, denoted
$\varepsilon^{2}$) is used as the loss function.

To make the boundary effects more visible, we divide each test patch into two
parts, interior and frame, with a frame width of eight pixels. The
reconstruction error is computed separately over these two parts. For our
patches of size $192\times 192$, $\varepsilon^{2}_{\text{INTER}}$ is computed
over the central part of size $176\times 176$ ($176=192-2\times 8$), and
$\varepsilon^{2}_{\text{FRAME}}$ over the cropped frame. The MSEs are
summarised in Table I, which show that our method (Diff) has not only achieved
the highest accuracy for both interior and frame but also maximally reduced
the error gap between interior and frame. We visualise the error maps over a
randomly picked region (near Caspian Sea), as shown in Fig. 3e$\sim$3m. These
error maps are obtained in three steps: patch reconstruction by the U-Net,
assembling the $4\times 4$ non-overlapping error maps, and applying the Farid
transform [38] to detect the horizontal and vertical edges. It is shown that
the boundary artefacts are visible in all the error maps except the one
delivered by Diff. Table I and Fig. 3 make it evident that Diff performs
significantly better than the other methods for this super-resolution task.

Figure 3: Super-resolution reconstruction of a world topographic map using a
U-Net. The top row shows the data at three scales, (a) for global, (b) zooming
into the region marked by the red-coloured cross in (a), and (c) zooming into
the red-coloured box in (b). (c) shows a pair of input and output of the
U-Net. (d) shows the reconstruction of (b) by model Zero; the reconstructions
by the other models look similar. (e)$\sim$(m) display the error maps over
(b), obtained by patch reconstruction, assembling the $4\times 4$ error maps
(non-overlapping), and the Farid transform (both horizontal and vertical) to
highlight the artefacts near the patch boundaries (the horizontal and vertical
stripes).

##  V Conclusions

We have presented a new padding-free method for size-keeping image
convolution. The central idea is to establish an equivalence between window-
wise convolution over a discrete image and pixel-wise differentiation over a
continuous representation of that image. Convolution within an incomplete
sliding window can then be achieved by differentiation at its centre, with the
continuous image parameterised from the nearest complete window. As such, our
method preserves the differential characteristics of kernels. Our final
formula is simple and computationally lightweight, available for both image
filtering and CNN-based machine learning. The preservation of the differential
operator at the boundary pixels makes our method more accurate for processing
images with smoother boundaries, such as mathematical or physical fields and
high resolution images. Our experiments have shown visible superiority of our
method on both image filtering and CNN-based computer vision tasks. We provide
an optimised implementation of our method, including both forward filtering
(to replace torch.nn.functional.conv2d) and a convolutional layer class (to
replace torch.nn.Conv2d), available open-source from https://github.com/stfc-
sciml/DifferentialConv2d (with all experiments included).

## Acknowledgements

This work is supported by the EPSRC grant, Blueprinting for AI for Science at
Exascale (BASE-II, EP/X019918/1), which is Phase II of the Benchmarking for AI
for Science at Exascale (BASE) grant.

## References

  * [1] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner, “Gradient-based learning applied to document recognition,” _Proceedings of the IEEE_ , vol. 86, no. 11, pp. 2278–2324, 1998.
  * [2] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [3] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” _Communications of the ACM_ , vol. 60, no. 6, pp. 84–90, 2017.
  * [4] G. Huang, Z. Liu, L. Van Der Maaten, and K. Q. Weinberger, “Densely connected convolutional networks,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2017, pp. 4700–4708.
  * [5] D. P. Kingma and M. Welling, “Auto-encoding variational bayes,” _arXiv preprint arXiv:1312.6114_ , 2013.
  * [6] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial networks,” _Communications of the ACM_ , vol. 63, no. 11, pp. 139–144, 2020.
  * [7] J. Ho, A. Jain, and P. Abbeel, “Denoising diffusion probabilistic models,” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 6840–6851, 2020.
  * [8] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole, “Score-based generative modeling through stochastic differential equations,” _arXiv preprint arXiv:2011.13456_ , 2020.
  * [9] S. Ioffe and C. Szegedy, “Batch normalization: Accelerating deep network training by reducing internal covariate shift,” in _International conference on machine learning_. pmlr, 2015, pp. 448–456.
  * [10] F. Chollet, “Xception: Deep learning with depthwise separable convolutions,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2017, pp. 1251–1258.
  * [11] C. Innamorati, T. Ritschel, T. Weyrich, and N. J. Mitra, “Learning on the edge: Explicit boundary handling in cnns,” _arXiv preprint arXiv:1805.03106_ , 2018.
  * [12] G. Liu, K. J. Shih, T.-C. Wang, F. A. Reda, K. Sapra, Z. Yu, A. Tao, and B. Catanzaro, “Partial convolution based padding,” _arXiv preprint arXiv:1811.11718_ , 2018.
  * [13] J. Liu, “Improved padding in cnns for quantitative susceptibility mapping,” _arXiv preprint arXiv:2106.15331_ , 2021.
  * [14] B. Alsallakh, N. Kokhlikyan, V. Miglani, J. Yuan, and O. Reblitz-Richardson, “Mind the pad–cnns can develop blind spots,” _arXiv preprint arXiv:2010.02178_ , 2020.
  * [15] M. A. Islam, M. Kowal, S. Jia, K. G. Derpanis, and N. D. Bruce, “Position, padding and predictions: A deeper look at position information in cnns,” _arXiv preprint arXiv:2101.12322_ , 2021.
  * [16] C. H. Lin, H.-Y. Lee, H.-Y. Tseng, M. Singh, and M.-H. Yang, “Unveiling the mask of position-information pattern through the mist of image features,” _arXiv preprint arXiv:2206.01202_ , 2022.
  * [17] H. Gao, L. Sun, and J.-X. Wang, “Phygeonet: Physics-informed geometry-adaptive convolutional neural networks for solving parameterized steady-state pdes on irregular domain,” _Journal of Computational Physics_ , vol. 428, p. 110079, 2021.
  * [18] ——, “Super-resolution and denoising of fluid flow using physics-informed convolutional neural networks without high-resolution labels,” _Physics of Fluids_ , vol. 33, no. 7, p. 073603, 2021.
  * [19] S. Shen, H. Lu, M. Sadoughi, C. Hu, V. Nemani, A. Thelen, K. Webster, M. Darr, J. Sidon, and S. Kenny, “A physics-informed deep learning approach for bearing fault detection,” _Engineering Applications of Artificial Intelligence_ , vol. 103, p. 104295, 2021.
  * [20] Z. Fang, “A high-efficient hybrid physics-informed neural networks based on convolutional neural network,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 33, no. 10, pp. 5514–5526, 2021.
  * [21] P. Ren, C. Rao, Y. Liu, J.-X. Wang, and H. Sun, “Phycrnet: Physics-informed convolutional-recurrent network for solving spatiotemporal pdes,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 389, p. 114399, 2022.
  * [22] A. Alguacil, W. G. Pinto, M. Bauerheim, M. C. Jacob, and S. Moreau, “Effects of boundary conditions in fully convolutional networks for learning spatio-temporal dynamics,” in _Machine Learning and Knowledge Discovery in Databases. Applied Data Science Track: European Conference, ECML PKDD 2021, Bilbao, Spain, September 13–17, 2021, Proceedings, Part V_. Springer, 2021, pp. 102–117.
  * [23] V. Gupta and N. Ramani, “A note on convolution and padding for two-dimensional data,” _Geophysical Prospecting_ , vol. 26, no. 1, pp. 214–217, 1978.
  * [24] F. Aghdasi and R. K. Ward, “Reduction of boundary artifacts in image restoration,” _IEEE Transactions on Image Processing_ , vol. 5, no. 4, pp. 611–618, 1996.
  * [25] O. S. Kayhan and J. C. v. Gemert, “On translation invariance in cnns: Convolutional layers can exploit absolute spatial location,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 14 274–14 285.
  * [26] J. C. Myburgh, C. Mouton, and M. H. Davel, “Tracking translation invariance in cnns,” in _Artificial Intelligence Research: First Southern African Conference for AI Research, SACAIR 2020, Muldersdrift, South Africa, February 22-26, 2021, Proceedings 1_. Springer, 2020, pp. 282–295.
  * [27] A.-D. Nguyen, S. Choi, W. Kim, S. Ahn, J. Kim, and S. Lee, “Distribution padding in convolutional neural networks,” in _2019 IEEE International Conference on Image Processing (ICIP)_. IEEE, 2019, pp. 4275–4279.
  * [28] N. Yang, L. Zhong, F. Huang, D. Yuan, and W. Bao, “Random padding data augmentation,” _arXiv preprint arXiv:2302.08682_ , 2023.
  * [29] Y.-H. Huang, M. Proesmans, and L. Van Gool, “Context-aware padding for semantic segmentation,” _arXiv preprint arXiv:2109.07854_ , 2021.
  * [30] S. Wu, G. Wang, P. Tang, F. Chen, and L. Shi, “Convolution with even-sized kernels and symmetric padding,” _Advances in Neural Information Processing Systems_ , vol. 32, 2019.
  * [31] H.-T. Cheng, C.-H. Chao, J.-D. Dong, H.-K. Wen, T.-L. Liu, and M. Sun, “Cube padding for weakly-supervised saliency prediction in 360 videos,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 1420–1429.
  * [32] M. Sokouti, A. Zakerolhosseini, and B. Sokouti, “Medical image encryption: an application for improved padding based ggh encryption algorithm,” _The open medical informatics journal_ , vol. 10, no. 1, 2016.
  * [33] G. Liu, A. Dundar, K. J. Shih, T.-C. Wang, F. A. Reda, K. Sapra, Z. Yu, X. Yang, A. Tao, and B. Catanzaro, “Partial convolution for padding, inpainting, and image synthesis,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2022.
  * [34] Z. Li, H. Zheng, N. Kovachki, D. Jin, H. Chen, B. Liu, K. Azizzadenesheli, and A. Anandkumar, “Physics-informed neural operator for learning partial differential equations.(2021),” _arXiv preprint arXiv:2111.03794_ , 2021\.
  * [35] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in _Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18_. Springer, 2015, pp. 234–241.
  * [36] L. Fei-Fei, R. Fergus, and P. Perona, “One-shot learning of object categories,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 28, no. 4, pp. 594–611, 2006.
  * [37] N. N. C. for Environmental Information, “Etopo 2022 15 arc-second global relief model,” 2023.
  * [38] H. Farid and E. P. Simoncelli, “Differentiation of discrete multidimensional signals,” _IEEE Transactions on image processing_ , vol. 13, no. 4, pp. 496–508, 2004.


* [32] J. Pamula, S. Jajodia, P. Ammann, and V. Swarup, “A weakest-adversary security metric for network configuration security analysis,” in _Proceedings of the 2nd ACM workshop on Quality of protection_ , 2006, pp. 31–38.
  * [33] S. Noel and S. Jajodia, “Metrics suite for network attack graph analytics,” in _Proceedings of the 9th Annual Cyber and Information Security Research Conference_ , 2014, pp. 5–8.
  * [34] L. Wang, A. Singhal, and S. Jajodia, “Measuring the overall security of network configurations using attack graphs,” in _IFIP Annual Conference on Data and Applications Security and Privacy_. Springer, 2007, pp. 98–112.
  * [35] L. Wang, T. Islam, T. Long, A. Singhal, and S. Jajodia, “An attack graph-based probabilistic security metric,” in _IFIP Annual Conference on Data and Applications Security and Privacy_. Springer, 2008, pp. 283–296.
  * [36] L. Munoz-González and E. C. Lupu, “Bayesian attack graphs for security risk assessment,” _ST-153/RWS-21 NATO Workshop on Cyber Resilience_ , pp. 64–77, 2017.


11institutetext: Institute of Physics and Astronomy, University of Potsdam,
14476 Potsdam, Germany 22institutetext: Dublin Institute for Advanced Studies,
Astronomy & Astrophysics Section, 31 Fitzwilliam Place, D02 XF86 Dublin 2,
Ireland 33institutetext: DESY, 15738 Zeuthen, Germany 44institutetext: Centre
for Space Research, North-West University, 2520 Potchefstroom, South Africa
55institutetext: Astronomical Observatory of Ivan Franko National University
of Lviv, Kyryla i Methodia 8, 79005 Lviv, Ukraine

# Morphology of  supernova remnants and their halos

R. Brose , Corresponding author, 1122<EMAIL_ADDRESS>M. Pohl 1133 I.
Sushch 4455

(Received ; accepted)

###### Abstract

Context. Supernova remnants are known to accelerate particles to relativistic
energies, on account of their non-thermal emission. The observational progress
from radio to gamma-ray observations reveals more and more morphological
features that need to be accounted for when modeling the emission from those
objects.

Aims. We use our time-dependent acceleration code RATPaC to study the
formation of extended gamma-ray halos around supernova remnants and the
morphological implications that arise when the high-energetic particles start
to escape from the remnant.

Methods. We performed spherically symmetric 1-D simulations in which we
simultaneously solve the transport equations for cosmic rays, magnetic
turbulence, and the hydrodynamical flow of the thermal plasma in a volume
large enough to keep all cosmic rays in the simulation. The transport
equations for cosmic-rays and magnetic turbulence are coupled via the cosmic-
ray gradient and the spatial diffusion coefficient of the cosmic rays, while
the cosmic-ray feedback onto the shock structure can be ignored. Our
simulations span 25,000 years, thus covering the free-expansion and the Sedov-
Taylor phase of the remnant’s evolution.

Results. We find strong difference in the morphology of the gamma-ray emission
from supernova remnants at later stages dependent on the emission process. At
early times both - the inverse-Compton and the Pion-decay morphology - are
shell-like. However, as soon as the maximum-energy of the freshly accelerated
particles starts to fall, the inverse-Compton morphology starts to become
center-filled whereas the Pion-decay morphology keeps its shell-like
structure. Escaping high-energy electrons start to form an emission halo
around the remnant at this time. There are good prospects for detecting this
spectrally hard emission with the future Cerenkov Telescope Array, as there
are for detecting variations of the gamma-ray spectral index across the
interior of the remnant. Further, we find a constantly decreasing non-thermal
X-ray flux that makes a detection of X-ray unlikely after the first few
thousand years of the remnants evolution. The radio flux is increasing
throughout the SNR’s lifetime and changes from a shell-like to a more center-
filled morphology later on.

###### Key Words.:

Supernova Remnants - Cosmic Rays - Magnetic Turbulence

##  1 Introduction

The detection of extended gamma-ray emission around two nearby pulsars -
Geminga and PSR B0656+14 (Abeysekara et al., 2017a, b) \- spawned an extensive
discussion of TeV-halos around pulsars (Giacinti et al., 2020).

Escaping high-energy electrons and positrons can be confined close the
accelerating pulsar-wind nebulae (PWN) by self-amplified magnetic turbulence
and enhance the TeV gamma-ray emission in the vicinity of the PWN (Evoli et
al., 2018). The resulting halos will have an energy-dependent morphology, with
a smaller extension towards the highest energies (Principe et al., 2020).

However, the interpretation of the observational data is not trivial. Simple
uniform-diffusion models can explain the observed morphology of the gamma-ray
emission on account of confining all electrons and positrons close to the
pulsar. A contribution to the positron excess observed in the AMS-data (Chang
et al., 2008; Adriani et al., 2009) would thus be ruled out. In reality, a
spatially non-uniform diffusion coefficient has to be expected for self-
amplified turbulence and could both explain the observed morphology and allow
for sufficient escape flux to support the local positron-flux (Profumo et al.,
2018).

A similar self-regulation of the diffusion coefficient by escaping particles
is known to exist around supernova remnants (SNRs) as well. There, the
scattering turbulence is created mainly by escaping hadrons, and the electrons
accelerated at the SNR blast-wave are trapped as a side effect (Ohira et al.,
2010; Nava et al., 2016).

So far, there has not been a direct measurement of  cosmic rays (CRs) escaping
from a SNR. A possible scenario would be CRs illuminating molecular clouds
close the SNRs where the enhanced target density boosts hadronic gamma-ray
emission (Yan et al., 2012). There is evidence of larger extension of the
gamma-ray emission around RXJ1713.7-3946 compared to the X-ray emission that
could indicate CR escape (H. E. S. S. Collaboration et al., 2018a). However,
there is indirect evidence that CR escape has to happen around SNRs. Recent
studies of particle acceleration in supernova remnants showed that typically
soft, broken power-law spectra of aged SNRs (Zeng et al., 2019) can be
produced by the escape of particles at the highest energy from the interior of
the SNR. Here, the most energetic particles escape once the SNR is not capable
of accelerating them any further, creating soft spectra inside the SNR where
most of the emission is produced (Brose et al., 2020; Celli et al., 2019).

These concepts potentially resolve the tension between  the $s=2$ spectra
predicted by shock-acceleration theory and the  typically soft emission
spectra  of $s\geq 2.7$ observed inside  evolved SNRs , especially in the
gamma-ray domain. However, the question whether SNRs are the sources of the
Galactic CRs is yet still unanswered. The earlier studies show that there is
no  contradiction between the observed soft emission spectra and the somewhat
harder total-production spectra  with $s\approx 2.2-2.4$ that are predicted
for the sources of Galactic CRs by Galactic propagation models (Trotta et al.,
2011). The accumulated proton spectrum is approximately the  same as the time
integral over the SNRs lifetime and represents the CR yield that gets finally
released into the sea of Galactic CRs once the SNR-shock faded.

The aim of this paper is to explore the observational signatures that can be
expected from escaping CRs around SNRs.

##  2 Basic equations and assumptions

The methodology is similar to that of Brose et al. (2020), and here we shall
only give a short summary. We combine a kinetic treatment of the CRs with a
thermal leakage injection model, a fully time-dependent treatment of the
magnetic turbulence, and a PLUTO-based simulation of the hydrodynamical flow
profiles.

###  2.1 Cosmic rays

We solve the kinetic equation for the differential number density of CRs, $N$,

$\displaystyle\frac{\partial N}{\partial t}=$ $\displaystyle\nabla(D_{r}\nabla
N-\mathbf{u}N)$ $\displaystyle-\frac{\partial}{\partial
p}\left((N\dot{p})-\frac{\nabla\cdot\mathbf{u}}{3}Np\right)+Q$ (1)

in the test-particle limit, where $D_{r}$ denotes the spatial diffusion
coefficient, u the advective velocity, $\dot{p}$ energy losses  (see section
2.1.2), and $Q$ the source of thermal particles (Skilling, 1975).

We rewrite equation (1) by transforming the radial coordinate, $r$, to a new
spatial coordinate that is co-moving with the SNR shock and provides a very
high numerical resolution at the shock front (Brose et al., 2020, and
references therin). The sea of Galactic CRs is neglected in this simulations.
As long as the fraction of injected particles is constant over time, the
emission from freshly injected particles is always dominant over the
contribution from background CRs.

####  2.1.1 Injection

We use a thermal leakage model (Blasi et al., 2005; Malkov, 1998) for the
injection of particles. Here, the efficiency of injection $\eta_{i}$ is given
by

$\displaystyle\eta_{i}=\frac{4}{3\sqrt{\pi}}(\sigma-1)\psi^{3}e^{-\psi^{2}}\text{
, }$ (2)

where $\sigma$ is the shock compression ratio, and $\psi$ is the multiple of
the thermal momentum, at which we inject particles.

Several authors noted that the bipolar morphology in the non-thermal emission
of SN 1006 can be attributed to effects of the shock-obliquity on the
injection or acceleration efficiency (Völk et al., 2003; Petruk et al., 2009;
Beshley & Petruk, 2012; Pais & Pfrommer, 2020). This notion seem to be
supported by the hybrid-simulations of Caprioli & Spitkovsky (2014), where the
acceleration at quasi-perpendicular shocks is strongly suppressed. Reville &
Bell (2013) on the other hand found a quasi-universal behaviour of shocks
irrespective of the magnetic-field orientation very far upstream of the shocks
by using a spherical-harmonics expansion of the CR Fokker-Planck equation.
This suggests that injection may only be weakly dependent on the shock
orientation. We have to ignore possible obliquity effects in our spherically
symmetric model and assume a quasi-parallel configuration across the shock-
surface.

This injection scenario is a simplification, in particular for electrons, for
which pre-acceleration to a few tens of MeV is required and established at the
shock (Matsumoto et al., 2017; Li et al., 2018; Bohdan et al., 2019). We are
interested in particles at energies well above $100$ MeV, and so the
particulars of that pre-acceleration can be ignored.

The multiple of the thermal momentum, $\psi$, determines the fraction of
thermal particles that get turned into CRs. However, $\psi$ and hence $\eta$
are only weakly constrained. Low values for $\psi$ are required in scenarios
featuring non-linear modifications of the shock structure by CR pressure. The
originally proposed lower boundary $\psi\approx 3.5$ ensures that thermal
particles  have a mean-free path larger than the shock-thickness and can
participate in the DSA-process. For this work we choose $\psi=4.2$, which
guarantees a CR pressure of less than $2.5\,$% of the shock ram pressure
during the entire simulation. Likewise, this value is close to the amount of
injected particles seen in SN 1006, where the likely leptonic origin of the
emission allows a reasonably sound estimate of $\eta_{i}$ (see section A for
details). We have chosen an equal amount of injected electrons and protons.

####  2.1.2 Inverse-Compton losses

In addition to the synchrotron losses  that were considered already in earlier
versions of RATPaC, electrons will also suffer losses by collisions with
photons from background photon-fields. Usually these inverse-Compton (IC)
losses can be neglected as synchrotron losses in amplified magnetic field will
dominate. However, high-energy electrons spend a large fraction of time
upstream of the shock, where IC and synchrotron losses are of similar
strength.

We used the approximations derived by Reimer et al. (2006) to account for IC
losses by collisions with the cosmic microwave background. We note that
additional photon fields might have to be taken into account when modeling
core-collapse supernova remnants where strong local infrared and optical
photon fields have to be expected or when specific Type-Ia SNRs are modeled
for which estimates of the background photon-fields are available based on
their location in the Galactic disk.

###  2.2 Magnetic turbulence

In parallel to the transport equation for CRs, we solve a transport equation
for the magnetic-turbulence spectrum, assuming Alfvén waves only. The temporal
and spatial evolution of the spectral energy-density per unit logarithmic
bandwidth, $E_{w}$, is described by

$\displaystyle\frac{\partial E_{w}}{\partial
t}+\cdot\nabla(\mathbf{u}E_{w})+k\frac{\partial}{\partial
k}\left({k^{2}}D_{k}\frac{\partial}{\partial k}\frac{E_{w}}{k^{3}}\right)=$
$\displaystyle=2(\Gamma_{g}-\Gamma_{d})E_{w}\text{ . }$ (3)

Here, $\mathbf{u}$ denotes the advection velocity, $k$ the wavenumber, $D_{k}$
the diffusion coefficient in wavenumber space, and $\Gamma_{g}$ and
$\Gamma_{d}$ the growth and damping terms, respectively (Brose et al., 2016).

We calculate the diffusion coefficient from $E_{w}$ using

$\displaystyle D_{r}$
$\displaystyle=\frac{4v}{3\pi}r_{g}\frac{U_{m}}{{E}_{w}}\text{ , }$ (4)

where $U_{m}$ denotes the energy density of the large-scale magnetic field,
$v$ is the particle velocity, and $r_{g}$ the gyro-radius of the particle.

As initial condition, we used a diffusion coefficient, and hence a turbulence
spectrum, as suggested by Galactic propagation modeling (Trotta et al., 2011),
but reduced by a factor ten on account of numerical constraints,

$\displaystyle D_{0}$
$\displaystyle=10^{28}\left(\frac{pc}{10\,\text{GeV}}\right)^{1/3}\left(\frac{B_{0}}{3\,\mu\text{G}}\right)^{-1/3}\text{
. }$ (5)

This choice is a factor of ten higher then in our previous works.

We use growth-rate based on the resonant streaming instability (Skilling,
1975; Bell, 1978),

$\displaystyle\Gamma_{g}$
$\displaystyle=A\cdot\frac{v_{\text{A}}p^{2}v}{3E_{\text{w}}}\left|\frac{\partial
N}{\partial r}\right|\text{ , }$ (6)

where $v_{\text{A}}$ is the Alfven-velocity. We introduced a linear scaling
factor, $A$, to artificially enhance the amplification. We used $A=10$
throughout this paper to mimic the more efficient amplification due to the
non-resonant streaming instability (Lucek & Bell, 2000; Bell, 2004). The
particulars of the non-resonant amplification are beyond the present
capabilities of our code since the back-reaction of CR streaming that
terminates the wave-growth, a modification of the bulk flow (Riquelme &
Spitkovsky, 2009; Niemiec et al., 2010; Kobzar et al., 2017), can not be
accounted for. Likewise difficult to handle, and in fact quite unclear, is the
scattering efficiency of the non-resonant modes. The value we choose for $A$
guarantees a cut-off energy in the gamma-ray spectrum between $1-10\,$TeV as
observed in young SNRs.

We calculate the total magnetic-field strength as

$\displaystyle B_{\text{tot}}$
$\displaystyle=\sqrt{B_{0}^{2}+B_{\text{Turb}}^{2}}\text{ , }$ (7)

where $B_{0}$ is the large-scale magnetic field. We solve the induction
equation to model the transport of the frozen-in large-scale magnetic field
(Telezhinsky et al., 2013). The far-upstream field is assumed to be uniform
with strength $5\mu$G  and the field is assumed to be fully turbulent,
resulting in a magnetic field-compression of
$B_{\text{d}}=\sqrt{11}B_{\text{u}}$.

Since we exceed the growth rate of the resonant streaming instability (Bell,
1978) by a factor of ten, the turbulent field is amplified to $\delta B\gg
B_{0}$ during the initial phases of SNR evolution. The peak amplitude of the
field is reached right at the shock. Downstream the field strength quickly
falls due to efficient cascading of turbulence. The resulting magnetic-field
profiles resemble the profiles suggested by Pohl et al. (2005). In the early
phases of SNR evolution peak amplitudes of $\approx 100\,\mu$G are reached,
whereas after $1,000\,$years one finds $\approx 40\,\mu$G - a value compatible
with the $30\,\mu$G estimated for the $1,000\,$year-old SNR SN1006 (Acero et
al., 2010).

The growth of the magnetic turbulence and hence the magnetic field is balanced
by cascading. This process is described as a diffusion process in wavenumber
space, and the diffusion-coefficient is given by (Zhou & Matthaeus, 1990;
Schlickeiser, 2002)

$\displaystyle D_{\text{k}}$
$\displaystyle=k^{3}v_{\text{A}}\sqrt{\frac{E_{\text{w}}}{2B_{0}^{2}}}\text{ .
}$ (8)

This phenomenological treatment will result in a Kolmogorov-like spectrum, if
cascading is dominant. Since $v_{\text{A}}\propto B_{\text{tot}}$, the
cascading rate will depend on the level of magnetic turbulence in two
different regimes

$\displaystyle D_{k}\propto\begin{cases}\sqrt{E_{\text{w}}}&\text{for
$E_{\text{w}}<B_{0}^{2}/8\pi$}\\\ E_{\text{w}}^{3/2}&\text{for
$E_{\text{w}}>B_{0}^{2}/8\pi$}\text{ . }\\\ \end{cases}$ (9)

Once the turbulent field dominates over the background field, the cascading
rate depends more sensitively on the energy density of magnetic turbulence.
Without requiring other damping mechanisms, the enhanced cascading efficiently
limits the maximum level of turbulence to a level commensurate with that
derived from SNR observations and PIC simulations (Vink, 2006; Riquelme &
Spitkovsky, 2009; Niemiec et al., 2010).

###  2.3 Thermal plasma

In the test-particle limit, the evolution of an SNR can be described with the
standard gas-dynamical equations:

$\displaystyle\frac{\partial}{\partial t}\left(\begin{array}[]{c}\rho\\\
\textbf{m}\\\
E\end{array}\right)+\nabla\left(\begin{array}[]{c}\rho\textbf{v}\\\
\textbf{mv}+P\textbf{I}\\\ (E+p)\textbf{v}\end{array}\right)^{T}$
$\displaystyle=\left(\begin{array}[]{c}0\\\ 0\\\ L\end{array}\right)$ (19)
$\displaystyle\frac{\rho\textbf{v}^{2}}{2}+\frac{P}{\gamma-1}$
$\displaystyle=E\text{,}$ (20)

where $\rho$ is the density of the thermal gas, v the plasma velocity,
$\textbf{m}=\textbf{v}\rho$ the momentum density, $P$ the thermal pressure of
the gas, $L$ the energy losses due to cooling, and $E$ the total energy
density of the ideal gas with $\gamma=5/3$. We solve this system of equations
under the assumption of spherical symmetry in 1-D using the PLUTO code
(Mignone et al., 2007). The non-equilibrium cooling function, $L$, is taken
from Sutherland & Dopita (1993).

In this work, we display results for type-Ia supernova explosions. We initiate
the simulations with exponential-ejecta profiles:

$\displaystyle\rho_{\mathrm{SN}}=$ $\displaystyle
A\exp(-v/v_{e})t_{i}^{-3,}\text{ and }v=r/t_{i},$ (21) with $\displaystyle
v_{e}=\left(\frac{E_{\mathrm{ex}}}{6M_{\mathrm{ej}}}\right)^{1/2,}\text{ and
}A=\frac{6^{3/2}}{8\pi}\frac{M_{\mathrm{ej}}^{5/2}}{E_{\mathrm{ex}}^{3/2}}$
(22)

as initial conditions (Dwarkadas & Chevalier, 1998). Here,
$t_{i}=2.5\,\mathrm{yrs}$ is the start time of our simulation,
$M_{\mathrm{ej}}=1.4M_{\mathrm{sol}}$ the ejecta mass,
$E_{\mathrm{ex}}=10^{51}\,$erg the explosion energy, and $r$ the spatial
coordinate. The density of the ambient medium was chosen to be
$0.4\,\text{cm}^{-3}$.

##  3 Results

We followed the evolution of the remnant for $25,000\,$years. The remnant
enters the Sedov-Taylor phase after $1,300\,$years111The use the initial
conditions Eq. (21-22) results in a continous, smooth variation of the
expansion parameter $m=V_{\textbf{sh}}/(R_{\text{sh}}\cdot t)$. After 1,300
yrs $m=0.47$, reasonabliy close to the asymptotic $m=0.4$ clasically expected
for the Sedov-phase. and would enter the post-adiabatic phase after
$35,000\,$years.

The total magnetic field reaches $90\,\mu$G after $300\,$yrs and drops to
$40\,\mu$G after $1000\,$yrs. After $10,000\,$yrs, the turbulent field is
weaker than the compressed large-scale field.

In the following, we first describe the morphology of the TeV-halo over the
lifetime of the remnant and then examine its detectability with today’s and
future gamma-ray experiments.

###  3.1 Halo-evolution

We calculated intensity maps for  inverse-Compton and Pion-decay (PD) emission
at three energies over the lifetime of the SNR. The results are presented in
Figure 1.

Figure 1: Maps of normalized surface brightness of a Type-Ia SNR after 300,
1000, 2,000 and 10,000 yrs (from top left to bottom right). Left hemispheres
are for IC emission and right hemispheres for PD emission.

Initially, both radiation mechanisms produce a shell-like morphology in all
energy bands. From roughly $1000\,$years on, the shell thickness of IC
emission exceeds that of PD radiation, because the latter is boosted by the
high gas density immediately downstream of the shock. Contrary, the IC
emission reflects only the distribution of electrons, and already after
$1000\,$years we notice IC emission outside the SNRs shell. The age of
$1000\,$years roughly corresponds to the time when the SNR luminosity peaks
regardless of the  gamma-ray emission mechanism, on account of the transition
to the Sedov-Taylor stage. For an ambient density of $0.4\,$cm-3 and an
average free-expansion velocity of $5000\,$km/s, five solar masses of material
will have been swept-up by the shock within the first 1000 years of
expansion222The initial shock-velocity might exceed $10,000\,$ km/s but it is
not constant even for the free-expansion phase.. The transition to the Sedov-
phase usually also marks the time of the highest maximum energy of particles
(Ptuskin & Zirakashvili, 2003). At later times, particles of the highest
energy start to escape from the remnant (Brose et al., 2020).

Even after $10,000\,$years, the PD morphology remains shell-like whereas that
of IC emission becomes center-filled. Even low-energy CRs propagated into the
center of the remnant, and the projection enhances the brightness towards the
center of the remnant in the IC channel. The low density of the thermal plasma
in the center keeps suppressing PD emission from the central region. As more
and more high-energy electrons escape with increasing age, an extensive halo
of $>100$-GeV electrons is formed around the remnant. There is also a faint
halo in the PD channel after $10,000\,$years, but it is much weaker than the
IC halo on account of the low gas density.

Synchrotron losses most strongly modify the distribution of high-energy
electrons. They cause the thinner IC shell at $1\,$TeV compared to lower
energies that is visible in the intensity maps. The magnetic field peaks at
the shock, and so synchrotron losses are strongest there. The same effect is
responsible for the relative smaller extend of the bright IC-region at
$2,000\,$years. In the unshocked ejecta the field strength is $0.2\,\mu$G,
making IC scattering the dominant energy-loss mechanism. Generally, the IC-
loss timescale exceeds that for synchrotron losses where the magnetic field
surpasses $5\,\mu$G.

Figure 2 shows the projected profiles of gamma-ray intensity at four stages of
the remnant evolution.

Figure 2: Projected intensity profiles normalized to the surface brightness at
the shock at four stages of SNR evolution.

It is clearly visible that the extension of the IC halo increases with energy,
in contrast to the case of PWNs. The size of the halo is not defined by a
balance between acceleration of the highest energetic electrons and their
synchrotron cooling, as all high-energetic particles have been accelerated at
early times. More important is that the high-energy particles experience the
largest diffusion coefficient and thus can fill a larger halo.

Figure 3: Normalized diffusion coefficient in the upstream region as function
of distance to the shock and the SNR age. The left panel is for a particle
energy of $100\,$GeV, the right panel for $3\,$TeV. The dashed lines
correspond to $1\,$%, $10\,$% and $50\,$% of $D_{\text{ISM}}$.

The morphology of the remnant also strongly depends on irregularities in the
ambient medium. Deviations from spherical symmetry may be caused if the
coherence length of the turbulent ambient magnetic field is larger than the
size of the remnant (Pais & Pfrommer, 2020). Similar effects have been
described earlier for the hadronic and leptonic emission morphology of
remnants expanding in a uniform ambient field (Petruk et al., 2009; Beshley &
Petruk, 2012). Our results qualitatively agree with these conclusions in
regions where the field is parallel to the shock or where the coherence-scale
of the ambient field is smaller than the size of the remnant. However,
accelerated diffusion into the center of the remnant with the onset of the
Sedov-stage produces a more center-filled morphology for the IC-emission than
can be obtained by models relying solely on advection.

A comparison of our predictions with measurements is difficult as there is no
firm detection of an emission halo around a SNR (see also section 3.3 for
details). However, extended gamma-ray emission around the two Pulsars Geminga
and PSR B0656+14 has been detected (Abeysekara et al., 2017a). The
measurements show a roughly exponential decrease of the surface-brightness
with increasing radial distance from the central pulsar. This behaviour is
also reproduced in simple models using one or two zones with spatially
constant diffusion coefficients around the pulsars (Di Mauro et al., 2020). We
obtain a similar behaviour in our halo-profiles (see Figure 2) for those
evolutionary stages during which CRs are still accelerated to the highest
energies. However, near the cut-off energy the profiles transition to a more
linear trend reflecting that those particles simply escape. This behaviour is
more pronounced for electrons, on account of the synchrotron cooling of the
high-energy electrons in the strong magnetic field inside the remnant.

###  3.2 Reduction of the diffusion coefficient

The escape of particles from the SNR and consequently the amplification of
turbulence change the diffusion coefficient in the vicinity of the SNR (Fujita
et al., 2010, 2011). The escape of high-energy particles from their
acceleration sites gained new attention with the detection of an extended halo
around the Pulsars Geminga and PSR B0656+14 (Abeysekara et al., 2017a).

Figure 3 shows the spatial variation of the diffusion coefficient relative to
the assumed diffusion coefficient in the ISM, which is set to 10% of the
conventional Galactic diffusion coefficient (Trotta et al., 2011). This choice
was made to keep the time-step in our simulations reasonably large, otherwise
the simulation of very old SNRs would be impossible.

For both electrons and ions the reduction of the diffusion coefficient is
strongest after $5,000\,$years ($10$%-level) and $2,000\,$years ($1$%-level),
around the time of, or shortly after, the transition from free expansion to
the Sedov-Taylor phase, when the maximum energy of the accelerated particles
starts to decrease. Then there are too few freshly accelerated particles at
high energy to sustain the level of turbulence, and consequently the diffusion
coefficient starts to increase. The main damping mechanism of turbulence is
cascading which, in contrast to other studies in the context of Pulsars (Evoli
et al., 2018), we treat not as a simple loss term but account for the energy
transfer to smaller scales. Consequently, the distributions in Figure 3 show
an energy-dependent time-evolution. The spatial extent is similar though,
because the suppression of the diffusion coefficient for $100\,$GeV-particles
is governed mainly by down-cascading from larger scales, i.e. by turbulence
driven by particles at higher energy. Turbulence driving by the low-energy
particles is only important very close to the shock.

To better understand the energy-dependent time evolution of upstream diffusion
after the end of the free-expansion phase, it is instructive to closely
inspect the cascading time (Schlickeiser, 2002),

$\displaystyle\tau(k)$
$\displaystyle=\frac{1}{v_{A}k}\sqrt{\frac{2U_{B}}{E_{w}}},$ (23)

which depends on the wave-number, $k$, the spectral-energy density of the
turbulent field, $E_{w}$, the Alfvén speed, $v_{A}$, and the energy density in
the large-scale magnetic field, $U_{B}$. At all wave-numbers, except at
$k_{\text{min}}$ resonant with the most energetic particles, turbulence energy
is cascaded from large scales (small $k$) to small scales (large $k$). So, in
a quasi-steady state, the level of turbulence at small scales is sustained by
a continuous influx from larger scales.

At an age of $4,000\,$ years, the cascading times for turbulence resonant with
$3\,$TeV ($100\,$GeV) particles $5\,$pc ahead of the shock is about
$500\,$years ($30\,$years). It takes $4,300$ years for the shock to arrive in
this region, suggesting that cascading is efficient. However, the full
spectral transport of turbulence is in many cases slower than the cascading
time suggests, implying that the variation in the level of turbulence arises
from a changing CR density gradient. As the freshly accelerated particles
become fewer and less energetic with time, particles escaping from deep
downstream become more important for turbulence driving.

During the initial growth of turbulence in the precursor, for the two particle
energies shown in Figure 3 the diffusion coefficient shows the same trend,
until the supply of freshly accelerated $3\,$TeV-particles is exhausted.
However, the precursor scale of $100\,$GeV-particles is significantly smaller
than for $3\,$TeV-particles, because the turbulence scattering the lower-
energy particles results from cascading of modes resonant with the more
energetic particles. Later, the escape of particles, that were trapped in the
interior of the remnant, becomes the main driver of turbulence, and in fact
the maximum extent of the turbulence precursor, and hence the reduced
diffusion coefficient, is reached well after the beginning of the Sedov-phase.

It is important to emphasize again that the diffusion coefficient for low-
energy particles is determined by the turbulence that the high-energy
particles provide. This has consequences beyond Supernova remnants, as for
example low-energy particles in the halos of PWNs will reside closer to their
acceleration sites than is suggested, if one treats cascading as a simple
damping mechanism (Evoli et al., 2018).

###  3.3 Detectability

Figure 4 shows the emission spectra expected from the SNR itself and the halo
at different times. We assumed a generic distance of $1\,$kpc to the remnant.
The contribution from the halo is shown twice, once as volume-integrated
emission from $r>r_{sh}$, and once only the component that in projection
appears to come from  beyond the projected radius of the SNR.  We model three-
dimensional objects, but we observe them in their 2D appearance in the sky.
Thus, part of the emission from the halo will in projection appear to come
from the SNR-interior, as the point of origin is in front of or behind the
remnant.

Figure 4: Comparison of emission spectra for IC (red) and PD (black) emission.
Emission from the SNR is filled-black and filled-red for PD and IC-emission
respectively. The halo-emission is filled-orange and filled-gray for IC and
PD-emission respectively at 300, 1000, 2,000 and 10,000yrs (from top left to
bottom right). The spectra including (and excluding) the projection effect
constitute the upper (lower) boundaries of the SNR emission and the lower
(upper) limit of the area for emission from the halo.

The evolution of the spectral energy distributions (SEDs) for PD-emission
clearly shows the mechanism for the spectral softening: as highly energetic
protons escape from the inside of the remnant, the low-energy particles remain
inside, leading to a soft spectrum (Brose et al., 2020). The high-energy
particles outside the SNR lack target material for the production of
significant emission, and hence hadronic halo emission is unlikely detectable
with current-generation gamma-ray observatories, even if the emission from the
remnant and the halo could be disentangled.

In case of the IC-emission, cooling becomes relevant after about
$1,000\,$years, and the spectra from the interior of the remnant are modified
by it. The halo emission from IC scattering is always much brighter than the
PD component. Overall, the flux from the halo is only 20% to 30% of that of
the SNR itself, making a detection possible only for the brightest known
Galactic SNRs. As low-energy particles reside closer to the shock, they
contribute little to the projected halo emission whose spectrum is
consequently harder than the overall emission from the remnant.

Figure 4 demonstrates that leptonic emission detected from an SNR by current-
generation gamma-ray observatories probably includes emission from the halo.
Most of the halo emission is produced close to the remnant, and in projection
it is very difficult to observationally separate it from emission from the
interior. It will probably require the tenfold higher sensitivity of the
future Cerenkov Telescope Array (CTA) to directly detect halo emission from
the brightest TeV-remnants. What may have been seen already in RXJ1713.7-3946
is that the leptonic halo causes larger extent of the very-high gamma-ray
emission compared to that of the X-ray emission.

A special case is the remnant SN 1006 that has a bipolar morphology of non-
thermal emission (Pye et al., 1981; Acero et al., 2010), possibly originating
from the variation of shock obliquity in a dominant large-scale magnetic field
(Völk et al., 2003). Regions with a quasi-perpendicular magnetic field are
associated with inefficient ion acceleration (Caprioli & Spitkovsky, 2014).
The apparent alignment of the bipolar morphology with the plane of the sky
means that the emission from SN 1006 is much less distorted by projection
effects than are more spherical remnants. If so, the contribution from the
halo to the high-energy gamma-ray emission could be a factor of a few higher
than that from other remnants of similar age. A direct measurement of the halo
emission might still be difficult due to the limited spatial resolution of
current-generation  Imaging Air Cerenkov telescopes (IACTs) and the low flux
from SN 1006. In any case, a large fraction of the most energetic electrons
resides in the upstream region (or halo). Due to the lower magnetic-field
strength in that region, they contribute relatively little to the X-ray
emission of SN 1006 (see also section 3.5) and a lot to the gamma-ray
emission. Figure 4 indicates that this could lead to a spectral discrepancy
between X-rays and high-energy gamma rays that is difficult to reproduce in a
one-zone model. In the past, such a finding for SN 1006 was interpreted as
effect of non-linear shock modification (Acero et al., 2010). We suggest that
this discrepancy is simply a natural result of halo emission.

###  3.4 Spectral index distribution

The angular resolution of gamma-ray observatories is key to a study of spatial
variations of the gamma-ray spectra (H. E. S. S. Collaboration et al., 2018a,
b; Humensky & VERITAS Collaboration, 2015). So far, no significant variation
across the remnant has been detected.

Figure 5: Maps of spectral-index variation in three energy bands at 300, 1000,
2,000 and 10,000yrs (from top left to bottom right), relative to the average
spectral index of the remnant (including the halo). Left hemispheres are
inverse-Compton emission and right hemispheres Pion-decay emission.

The spectrum of emission from outside the remnant is typically harder than
that from the interior. Figure 5 illustrates spectral index variations across
the remnant at different ages. For emission from the interior of the remnant,
the deviations are moderate . Depending on the energy-band, $|\Delta s|$ can
reach values of up to $0.8$, and in general the emission from the interior of
the remnant is softer than the average spectral index from the remnant.
However, $|\Delta s|\lesssim 0.2$ is usually fulfilled for all but the highest
energies of IC-emission. There, the mean spectral index is close to that seen
upstream of the shock due to the large contribution of upstream-electrons to
the gamma-ray flux. Still, the spectral index is quite uniform across the
interior of the remnant, whereas the deviations of the spectral index of the
halo emission are typically larger.

The detection of hadronic halo emission is unlikely given the current
observational sensitivities. In the H.E.S.S. data for RXJ1713.7-3946 the
$1$-$\sigma$ statistical uncertainty in the spatially-resolved spectral index
maps is $\Delta s\approx\pm 0.1$ (H. E. S. S. Collaboration et al., 2018a),
likewise $\Delta s\approx\pm 0.15$ for Vela Junior (H. E. S. S. Collaboration
et al., 2018b). The systematic and the statistical uncertainties are
comparable. In the VERITAS results for IC443, the spectral uncertainty ranges
between $\Delta s=\pm 0.11$ and $\Delta s=\pm 0.42$ (Humensky & VERITAS
Collaboration, 2015). Considering the limited spatial resolution and the poor
statistics in some regions, a significant detection of spectral-index
deviations is unlikely with current-generation IACTs.

The higher sensitivity of CTA should reduce the uncertainty, $\Delta s$, by
roughly a factor of three, which should allow the detection of spectral-index
variation with $3\sigma$ significance. We note that deviations from a
spherical symmetry or asymmetric environments, such as density gradients in
the ISM or fast motion of the progenitor (Meyer et al., 2020, 2021), could
lead to additional significant variations in the measured spectral index
across the remnant.

###  3.5 Synchrotron emission

We evaluated the non-thermal synchrotron emission from the SNR over its
lifetime. Figure 6 shows the time-evolution of the total synchrotron flux, as
well as the flux from the downstream region only.

Figure 6: Synchrotron flux from the remnant at various times. The upper
boundaries of the filled areas mark the total emission from the remnant while
the lower boundaries indicate emission from downstream only. The green band
indicates the $1-10\,$GHz range.

Contrary to the IC-emission, there is little difference between the total and
the downstream-only flux, except for very old remnants, because a strong
magnetic field boosts the emission from the downstream. At later times, most
of the high-energy electrons were able to escape the remnant and the magnetic
field is not amplified any more, the halo emission starts to dominate.

The flux in the X-ray band decreases with time, because the magnetic-field
strength decreases and cooling takes its toll in the electron spectra.
Consequently, the remnant is brightest and best detectable during the first
$2,000\,$yrs after the supernova - matching the census of Galactic SNRs (Vink,
2012).

In contrast, the radio flux increases throughout  the simulations. Energy
losses are not important, and the steady accumulation of GeV-scale particles
compensates for the weakening magnetic field. Interestingly, the escape of
electrons affects the spectral index in the radio band. Between 1 GHz and 30
GHz (Roughly $10^{-8}\,$keV to $10^{-7}\,$keV) the radio spectra are fairly
soft after $10,000\,$yrs, with a spectral index $\alpha\approx 0.75$
($S_{\nu}\propto\nu^{-\alpha}$).

Figure 7 shows the emission morphology at radio and X-ray energies for four
stages of SNR evolution.

Figure 7: Maps of normalized synchrotron surface brightness of a Type-Ia SNR
after 300, 1000, 2,000 and 10,000 yrs (from top left to bottom right).

Throughout the SNR’s lifetime the magnetic-field strength peaks right
downstream of the forward shock. Initially, when the turbulent magnetic field
exceeds the compressed large-scale field, cascading causes a decline toward
the interior. Later, the trend persists due to the flow structure in the
downstream of the shock. Consequently, the synchrotron emissivity peaks
immediately downstream to the shock-location.

The X-ray morphology consists of a thin shell throughout the entire lifetime
even though a part of the soft X-ray emission originates from the halo after
$2000\,$years and might be detectable for bright sources. The cutoff of the
synchrotron spectrum is in or slightly below the X-ray band, and so the X-ray
output is very sensitive to the structure of the magnetic field (Pohl et al.,
2005). A slight shift in the characteristic synchrotron frequency due to a a
spatial variation of the magnetic field can impose a strong variation in the
synchrotron flux.

The radio emission shows the same shell-like structure. However, the spatial
variation of the magnetic field has a moderate impact on the synchrotron
emissivity, typically $\propto B^{1.5}$. The radio shells appear thicker as a
consequence. Moreover, diffusion of electron towards the center of the remnant
partially compensates the weak magnetic field in the interior. Thus, the
remnant appears somewhat center-filled in the radio band at later stages.
Additionally, 4% of the radio-emission at 10,000 years originates from the
halo. The magnetic field is dominated by the large-scale field at that time,
which boosts the radio-emission from a downstream electron by a factor of $6$
compared to that of an upstream electron. This means that $20$% of the low-
energy electrons reside upstream of the shock and are escaping the remnant
already.

To explain the center-filled radio morphology of the remnant, Sushch & Hnatyk
(2014) argued that additional, inward-moving shocks are needed to almost
uniformly fill the interior of the Vela SNR with radio-emitting electrons. Our
findings suggest that the same can be achieved by diffusive transport of
electrons, in particular when in later stages of evolution the diffusion
coefficient inside and outside of the remnant increases. However, Vela is
likely the result of a core-collapse SN and thus additional factors need to be
taken into account. Core-collapse SNRs initially evolve in the freely
expanding wind of the progenitor star, in which the magnetic field will follow
a $B\propto 1/r$ dependence. Beyond the wind-termination shock, the magnetic-
field strength is almost constant out to the contact discontinuity (Sushch et
al. (2018) and Sushch et al. (in preparation)), which strongly affects the
X-ray and radio maps and likely enhances the trend towards a center-filled
radio morphology.

##  4 Conclusions

We performed numerical simulations of particle acceleration in SNRs, solving
time-dependent transport equations of CRs and magnetic turbulence in the test-
particle limit alongside the standard gas-dynamical equations for Type-Ia
SNRs. We derived the CR diffusion coefficient from the spectrum of magnetic
turbulence, that evolves through driving by the resonant streaming instability
as well as cascading and wave damping.

The gamma-ray morphology of SNRs depends on the radiation mechanism. Leptonic
gamma-ray emission tends to produce a center-filled appearance with increasing
remnant age, whereas hadronic gamma-ray emitters retain a shell-like
morphology, on account of the distribution of target material for p-p
interactions.

We showed that a gamma-ray halo has to be expected for leptonic emitters,
despite the projection effects that place part of the halo emission inside the
projected SNR radius. Whereas a firm detection of halo emission is in sight
for the leptonic case, the halo component of hadronic emission likely remains
undetectable, unless additional target material is present in the vicinity of
the SNR.

Our simulations show that the most energetic particles tend to reside outside
the SNR. The relatively weak magnetic field in the upstream region leads to a
low intensity of X-ray emission from these high-energetic electrons, making
the X-ray spectrum appear softer than expected for a uniform magnetic field.
It is likely that the observed spectral between X-ray and gamma-ray emission
from SN 1006 is based on this effect instead of non-linear shock-modification.

Upstream of the forward shock, the diffusion coefficient is reduced, confining
the accelerated particles close to the SNR. There is a strong spatial and
temporal evolution of the diffusion coefficient in the entire precursor. The
diffusion coefficient of low-energy particles is determined by cascading of
turbulence downward from larger scales. In contrast to PWNs, for which neither
the abundance of relativistic ions nor the time profile of particle
acceleration are known, for SNRs we can self-consistently determine the
injection spectra of CRs that will create the region of reduced diffusion
around SNRs.

An investigation of the distribution of gamma-ray spectral indices across the
interior of SNRs showed that no deviation beyond a $|\delta s|\approx 0.2$ has
to be expected for the known and resolvable Galactic SNRs. Detecting this
small deviations is beyond the capabilities of current-generation IACTs but
might be possible with CTA. Strong deviations from spherical symmetry can
further enhance the detectability.

Initially the remnant emits non-thermal X-rays above $0.1\,$keV but cooling
and a decreasing magnetic field reduce the X-ray luminosity and make a
detection unlikely after a few thousand years. The radio flux from the remnant
increases throughout its lifetime, and the radio morphology evolves from
shell-like to a more center-filled appearance.

## Acknowledgements

Robert Brose acknowledges funding from an Irish Research Council Starting
Laureate Award (IRCLA/2017/83).

## References

  * Abeysekara et al. (2017a) Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2017a, Science, 358, 911
  * Abeysekara et al. (2017b) Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2017b, ApJ, 843, 40
  * Acero et al. (2010) Acero, F., Aharonian, F., Akhperjanian, A. G., et al. 2010, A&A, 516, A62
  * Adriani et al. (2009) Adriani, O., Barbarino, G. C., Bazilevskaya, G. A., et al. 2009, Nature, 458, 607
  * Bell (1978) Bell, A. R. 1978, MNRAS, 182, 147
  * Bell (2004) Bell, A. R. 2004, MNRAS, 353, 550
  * Beshley & Petruk (2012) Beshley, V. & Petruk, O. 2012, MNRAS, 419, 1421
  * Blasi et al. (2005) Blasi, P., Gabici, S., & Vannoni, G. 2005, MNRAS, 361, 907
  * Bohdan et al. (2019) Bohdan, A., Niemiec, J., Pohl, M., et al. 2019, ApJ, 878, 5
  * Brose et al. (2020) Brose, R., Pohl, M., Sushch, I., Petruk, O., & Kuzyo, T. 2020, A&A, 634, A59
  * Brose et al. (2016) Brose, R., Telezhinsky, I., & Pohl, M. 2016, A&A, 593, A20
  * Caprioli & Spitkovsky (2014) Caprioli, D. & Spitkovsky, A. 2014, ApJ, 783, 91
  * Celli et al. (2019) Celli, S., Morlino, G., Gabici, S., & Aharonian, F. A. 2019, MNRAS, 490, 4317
  * Chang et al. (2008) Chang, J., Adams, J. H., Ahn, H. S., et al. 2008, Nature, 456, 362
  * Di Mauro et al. (2020) Di Mauro, M., Manconi, S., & Donato, F. 2020, Phys. Rev. D, 101, 103035
  * Dwarkadas & Chevalier (1998) Dwarkadas, V. V. & Chevalier, R. A. 1998, ApJ, 497, 807
  * Evoli et al. (2018) Evoli, C., Linden, T., & Morlino, G. 2018, Phys. Rev. D, 98, 063017
  * Fujita et al. (2010) Fujita, Y., Ohira, Y., & Takahara, F. 2010, ApJ, 712, L153
  * Fujita et al. (2011) Fujita, Y., Takahara, F., Ohira, Y., & Iwasaki, K. 2011, MNRAS, 415, 3434
  * Giacinti et al. (2020) Giacinti, G., Mitchell, A. M. W., López-Coto, R., et al. 2020, A&A, 636, A113
  * H. E. S. S. Collaboration et al. (2018a) H. E. S. S. Collaboration, Abdalla, H., Abramowski, A., et al. 2018a, A&A, 612, A6
  * H. E. S. S. Collaboration et al. (2018b) H. E. S. S. Collaboration, Abdalla, H., Abramowski, A., et al. 2018b, A&A, 612, A7
  * Humensky & VERITAS Collaboration (2015) Humensky, B. & VERITAS Collaboration. 2015, in International Cosmic Ray Conference, Vol. 34, 34th International Cosmic Ray Conference (ICRC2015), 875
  * Kobzar et al. (2017) Kobzar, O., Niemiec, J., Pohl, M., & Bohdan, A. 2017, MNRAS, 469, 4985
  * Li et al. (2018) Li, R., Zhou, C. T., Huang, T. W., et al. 2018, Physics of Plasmas, 25, 082103
  * Lucek & Bell (2000) Lucek, S. G. & Bell, A. R. 2000, MNRAS, 314, 65
  * Malkov (1998) Malkov, M. A. 1998, Phys. Rev. E, 58, 4911
  * Matsumoto et al. (2017) Matsumoto, Y., Amano, T., Kato, T. N., & Hoshino, M. 2017, Phys. Rev. Lett.
  * Meyer et al. (2020) Meyer, D. M. A., Petrov, M., & Pohl, M. 2020, MNRAS, 493, 3548
  * Meyer et al. (2021) Meyer, D. M. A., Pohl, M., Petrov, M., & Oskinova, L. 2021, MNRAS, 502, 5340
  * Mignone et al. (2007) Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228
  * Nava et al. (2016) Nava, L., Gabici, S., Marcowith, A., Morlino, G., & Ptuskin, V. S. 2016, MNRAS, 461, 3552
  * Niemiec et al. (2010) Niemiec, J., Pohl, M., Bret, A., & Stroman, T. 2010, ApJ, 709, 1148
  * Ohira et al. (2010) Ohira, Y., Murase, K., & Yamazaki, R. 2010, A&A, 513, A17
  * Pais & Pfrommer (2020) Pais, M. & Pfrommer, C. 2020, MNRAS, 498, 5557
  * Petruk et al. (2009) Petruk, O., Beshley, V., Bocchino, F., & Orlando, S. 2009, MNRAS, 395, 1467
  * Pohl et al. (2005) Pohl, M., Yan, H., & Lazarian, A. 2005, ApJ, 626, L101
  * Principe et al. (2020) Principe, G., Mitchell, A. M. W., Caroff, S., et al. 2020, arXiv e-prints, arXiv:2006.11177
  * Profumo et al. (2018) Profumo, S., Reynoso-Cordova, J., Kaaz, N., & Silverman, M. 2018, Phys. Rev. D, 97, 123008
  * Ptuskin & Zirakashvili (2003) Ptuskin, V. S. & Zirakashvili, V. N. 2003, A&A, 403, 1
  * Pye et al. (1981) Pye, J. P., Punds, K. A., Rolf, D. P., et al. 1981, MNRAS, 194, 569
  * Reimer et al. (2006) Reimer, A., Pohl, M., & Reimer, O. 2006, ApJ, 644, 1118
  * Reville & Bell (2013) Reville, B. & Bell, A. R. 2013, MNRAS, 430, 2873
  * Riquelme & Spitkovsky (2009) Riquelme, M. A. & Spitkovsky, A. 2009, ApJ, 694, 626
  * Schlickeiser (2002) Schlickeiser, R. 2002, Cosmic Ray Astrophysics
  * Skilling (1975) Skilling, J. 1975, MNRAS, 172, 557
  * Sushch et al. (2018) Sushch, I., Brose, R., & Pohl, M. 2018, A&A, 618, A155
  * Sushch & Hnatyk (2014) Sushch, I. & Hnatyk, B. 2014, A&A, 561, A139
  * Sutherland & Dopita (1993) Sutherland, R. S. & Dopita, M. A. 1993, ApJS, 88, 253
  * Telezhinsky et al. (2013) Telezhinsky, I., Dwarkadas, V. V., & Pohl, M. 2013, A&A, 552, A102
  * Trotta et al. (2011) Trotta, R., Jóhannesson, G., Moskalenko, I. V., et al. 2011, ApJ, 729, 106
  * Vink (2006) Vink, J. 2006, in ESA Special Publication, Vol. 604, The X-ray Universe 2005, ed. A. Wilson, 319
  * Vink (2012) Vink, J. 2012, A&A Rev., 20, 49
  * Völk et al. (2003) Völk, H. J., Berezhko, E. G., & Ksenofontov, L. T. 2003, A&A, 409, 563
  * Yan et al. (2012) Yan, H., Lazarian, A., & Schlickeiser, R. 2012, ApJ, 745, 140
  * Zeng et al. (2019) Zeng, H., Xin, Y., & Liu, S. 2019, ApJ, 874, 50
  * Zhou & Matthaeus (1990) Zhou, Y. & Matthaeus, W. H. 1990, J. Geophys. Res., 95, 14881

##  Appendix A Electron injection derived from SN1006

The injection fraction, $\eta_{i}$, is a crucial parameter for our
simulations. The normalization of the proton spectrum determines the gradient
of the CR distribution in the upstream region and hence the growth rate of
magnetic turbulence. At late times, $\eta_{i}$ can be understood as the number
ratio of particles injected into shock acceleration and the total number of
thermal particles in the system,

$\displaystyle\eta_{i}$
$\displaystyle=\frac{\int_{E_{i}}^{\infty}N(E)\text{d}E}{\int\rho_{0}/m\text{d}V}\text{
, }$ (24)

where $E_{i}$ is the injection energy, $N(E)$ the differential number density
of CRs, and $\rho_{0}/m$ the number density of the ambient medium. If
$\eta_{i}$ is known, equation (2) can be used to determine $\psi$ for the
setup of the simulation.

$\eta_{i}$ can also be determined from SNR observations. We pointed out that
the luminosity-age distribution of Galactic SNRs can be related to the total
number of accelerated electrons in these remnants, if one considers the
density profile for the Sedov solution. Then, the gamma-ray emission above
$1\,$TeV is dominated by the IC channel (Brose et al. 2020).

If the ambient density is very low, as it is for SN 1006, $\eta_{i}$ can also
be inferred directly, because it is likely that all of the high-energy photons
are produced via IC emission. The analysis of the H.E.S.S.-collaboration
revealed a total energy in electrons of $W_{e}=3.3\cdot 10^{47}\,$erg, an
electron spectral-index of $s_{e}=2.1$, and a electron cutoff-energy of
$E_{c}=10\,$TeV (Acero et al. 2010). These values can now be used to determine
the total number of electrons present in SN 1006. We assume for simplicity
that energetic electrons are present with constant spectral index between
$E_{0}=10\,$MeV and $E_{c}=10\,$TeV.

$\displaystyle W_{e}$ $\displaystyle=\int_{E_{0}}^{E_{c}}E\cdot N(E)\text{d}E$
(25) $\displaystyle N(E)$ $\displaystyle=N_{0}\cdot E^{-s}\text{ for
$E_{0}<E<E_{c}$}$ (26)
$\displaystyle\frac{W_{e}}{\int_{E_{0}}^{E_{c}}N(E)\text{d}E}$
$\displaystyle=\frac{E_{c}^{2-s}-E_{0}^{2-s}}{E_{c}^{1-s}-E_{0}^{1-s}}\frac{1-s}{2-s}=\alpha$
(27) $\displaystyle\int_{E_{0}}^{E_{c}}N(E)\text{d}E$
$\displaystyle=\frac{W_{e}}{\alpha}=6.4\cdot 10^{51}$ (28)

The total number of thermal particles can be calculated by integrating the
density of the ambient medium over the volume filled by the SNR. SN 1006 has
an reported radius of $R=9.6\,$pc in a medium with $\rho_{0}/m=0.085\,$cm-3.

$\displaystyle\int\frac{\rho_{0}}{m}\text{d}V$ $\displaystyle=\frac{4}{3}\pi
R^{3}\frac{\rho_{0}}{m}=9.26\cdot 10^{57}\text{ . }$ (29)

This gives $\eta_{i}=6.9\cdot 10^{-7}$ for electrons. A correction is in
order, because not the entire shock surface of SN 1006 is efficiently
accelerating particles. Estimating that only $\leavevmode\nobreak\ 20\,$% of
SN 1006’s surface are contributing to the high-energy emission, we find
$\eta_{i}=3.5\cdot 10^{-6}$, which corresponds to $\psi=4.205$. This is very
close to the value $\psi=4.2$ that we used in our simulations.

Experimental data suggests that for particles arriving at earth the electron-
to-proton ratio is $K_{ep}=0.01$. We inject $1.2$ times more electrons than
protons to account for a composition with $10\,$% helium. Since electrons and
protons are injected at the same $\psi$ but at different energies, this
results in $K_{ep}\approx 0.01$.