# End-to-End Segmentation via Patch-wise Polygons Prediction

Tal Shaharabany
Tel-Aviv University
<EMAIL_ADDRESS>Lior Wolf
Tel-Aviv University
<EMAIL_ADDRESS>
###### Abstract

The leading segmentation methods represent the output map as a pixel grid. We
study an alternative representation in which the object edges are modeled, per
image patch, as a polygon with $k$ vertices that is coupled with per-patch
label probabilities. The vertices are optimized by employing a differentiable
neural renderer to create a raster image. The delineated region is then
compared with the ground truth segmentation. Our method obtains multiple
state-of-the-art results: 76.26% mIoU on the Cityscapes validation, 90.92% IoU
on the Vaihingen building segmentation benchmark, 66.82% IoU for the MoNU
microscopy dataset, and 90.91% for the bird benchmark CUB. Our code for
training and reproducing these results is attached as supplementary.

##  1 Introduction

Accurate segmentation in dense environments is a challenging computer vision
task. Many applications, such as autonomous driving [33], drone navigation
[25], and human-robot interaction [24], rely on segmentation for performing
scene understanding. Therefore, the accuracy of the segmentation method should
be high as possible.

In recent years, fully convolutional networks based on encoder-decoder
architectures [21], in which the encoder is pre-trained on ImageNet, have
become the standard tool in the field. Several techniques were developed to
better leverage the capacity of the architectures. For example, skip
connections between the encoder and the decoder were added to overcome the
loss of spatial information and to better propagate the training signal
through the network [27]. What all these methods have in common is the
representation of the output segmentation map as a binary multi-channel image,
where the number of channels is the number of classes $C$. In our work, we add
a second type of output representation, in which, for every patch in the image
and for every class, a polygon comprised of $k$ points represents the binary
mask.

Alternative image segmentation architectures are presented in Fig. 1. The most
common is the encoder-decoder one. The active contour methods, similar to our
work, employ a polygon representation. Such methods, however, are iterative,
and produce a single polygon, lacking the flexibility to cope with occlusions.
Transformer based network employ patch-wise self-attention maps to calculate a
spatial representation that attends globally to multiple parts of the image.
The output mask is then generated by an MLP decoder. Our method differs from
these approaches and produces multiple polygons, which are rendered by a
neural renderer to produce the final segmentation mask.

For optimization purposes, we use a neural renderer to translate the polygons
to a binary raster-graphics mask. The neural renderer provides a gradient for
the polygon vertices and enables us to maximize the overlap between the output
masks and the ground truth. Since there are polygons for every patch and every
class, whether the object appears there or not, we multiply the binary mask
with a low-resolution binary map that indicates the existence of each object.

The method achieves state-of-the-art results on four segmentation benchmarks,
Cityscapes, CUB-200-2011, Nucleus segmentation challenges (MoNuSeg) and on the
Vaihingen building segmentation benchmark.

|  |
---|---|---
(a)
(b) | (c) | (d)
Figure 1: Comparison of the proposed method to the common segmentation
methods. (a) The prevalent encoder-decoder architecture predicts a pixel-wise
segmentation mask. (b) Deep active contour approaches, in which the trained
network predicts the contour of the object as a single polygon. The vertices
of the polygon are updated iteratively. (c) Transformer-based segmentation
models first divide the input image into per-patch tokens, each associated
with a positional encoding. The transformer algorithm then adds multiple
layers of interactions between the tokens. An MLP model then generates a local
mask per each token. (d) In our solution, the network predicts a single
polygon for every patch in the image, in addition to low-resolution
probability maps (not shown). Each polygon is rendered into a local mask
segmentation mask, and all the masks are aggregated to obtain the final
segmentation mask.

##  2 Related work

Most modern approaches to image segmentation perform image-to-image mapping
with a Fully Convolutional Network [21]. Many encoder-decoder architectures
were proposed, e.g., [12, 32, 30, 14], most of which employ U-Net based skip
connections between parallel blocks of the encoder and the decoder [27].
Improved generalization is obtained by reducing the gap between the encoder
and decoder semantic maps [41]. Another contribution adds residual connections
for each encoder and decoder block and divides the image into patches with a
calculated weighted map for each patch as input to the model. [37]. Wang et
al. add attention maps to each feature map in the encoder-decoder block in
order to enlarge the representation between far areas in the image with good
computational efficiency [36]. This is extended by adding scalar gates to the
attention layers [34]. The same contribution also adds a Local-Global training
strategy (LoGo), which further improves performance. Dilated convolutions [38]
are used by many methods to preserve the spatial size of the feature map and
to enlarge the receptive field [5, 40, 39]. The receptive field can also be
enlarged by using larger kernels [26] and enriched with contextual
information, by using kernels of multiple scales [6]. Hardnet [4] employs a
Harmonic Densely Connected Network that is shown to be highly effective in
many real-time tasks, such as classification and object detection.

Polygon-based segmentation  Active contour, also known as snakes [15], was
used on many instance segmentation applications over the years. Such methods
iteratively minimize an energy function that includes both edge-seeking terms
and a shape prior until it stops [17, 8, 3]. The first deep learning-based
snakes predict a correction vector for each point, by considering the small
enclosing patch, in order to converge to the object outline [29]. Acuna et
al.[1] propose automated and semi-automated methods for object annotating and
segmentation based on Graph neural networks and reinforcement learning. Ling
et al.[20] used a Graph Convolutional Network (GCN) in order to fit the
polygon’s positions to the outline. A different approach predicts a shift map
and iteratively changes the positions of the polygon vertices, according to
this map [10]. In order to obtain a gradient signal, the authors use a neural
renderer that transforms the 2D polygon into a 2D segmentation map. Our method
also employs a neural renderer. However, we perform local predictions per
patch and recover the vertex locations of the polygon in each patch non-
iteratively. This allows us to perform finer local predictions, and overcome
discontinuities caused by occlusions or by the existence of multiple objects.

A recent contribution [19], refines the output instance mask of pretrained
Mask-RCNN algorithm by learning to correct the coordinates of the extracted
polygons. This algorithm obtains good performances for the cityscapes instance
segmentation benchmark but uses different settings than the rest of the
literature. For example, the entire image is used instead of a specified
bounding box, and the bounding box expansion differs from what is
conventionally used in published work. Since no code was published, we cannot
perform a direct comparison. From the algorithmic perspective, the methods
greatly differ in the type of polygons used (boundary vs. local), the type of
optimization (Chamfer Distance loss [13] on polygon coordinates vs. using a
neural renderer), and almost every other aspect.

Neural Renderers  The ground truth segmentation is given as an image, and in
order to compare the obtained solution with it and compute a loss term, it is
required to have a differentiable way to convert the vertex representation to
a rasterized form. This conversion task, when implemented as a neural network,
is referred to as “neural renderers”. On the backward pass, fully differential
neural renderers convert image-based loss to gradients on the vertices’
coordinates. Loper et al. [22] approximate this gradient-based on image
derivatives, while Kato et al. [16] approximates the propagated error,
according to the overlap between the rasterization output and the ground
truth.

|  |
---|---|---
(a) | (b) | (c)
Figure 2: Example of the process of combining the polygon data with the grid-
based predictions, in order to obtain a refined segmentation map. (a) The low-
resolution segmentation map $M$ is shown for a single object (zoom-in from the
Cityscapes dataset). White pixels denote in-object, blacks are for the
background. (b) The corresponding polygons $R_{p}$. Despite using $k=5$, many
polygons match the patch boundaries. Different polygons are drawn in different
colors and are slightly scaled to visualize them separately. (c) The combined
segmentation map $M_{o}$, in the final resolution. The polygons on the top
right part of $R_{p}$ are gated out by $M$.

##  3 Methods

The input to our method is an RGB image $I$, of dimensions ${3\times H\times
W}$. The image is divided into patches of size $s\times s$ and a polygon with
$k$ vertices is predicted for each patch. In addition, each patch has an
associate scalar that determines the probability of the object being present
in the patch.

This mask representation is illustrated in Fig. 2. The low-res probability map
in panel (a) assigns one value per grid cell, where we denote higher
probabilities by a brighter value. The local polygons in panel (b) determine a
local segmentation mask in each grid cell. When the grid cell is completely
enclosed by the segmented object, the polygon covers the entire cell. At
object boundaries, it takes the form of the local boundary. When out of the
object, it can degenerate to a zero area polygon or an arbitrary one. The
final segmentation map, depicted in Fig. 2(c) is obtained by multiplying the
low-resolution probability map with the local patch rastered for each polygon.

The learned network we employ is given as $f=f_{2}\circ f_{1}$, where
$f_{1}:\mathbb{R}^{C\times{H}\times{W}}\rightarrow\mathbb{R}^{C\times\frac{H}{32}\times\frac{W}{32}}$
is the backbone (encoder) network, and $f_{2}$ is the decoder networks. The
output space of $f_{2}$ is a tensor that can be viewed as a field with
$\frac{H}{s}\times\frac{W}{s}$ vectors $\mathbb{R}^{2k+1}$, each representing
a single polygon that is associated with a single image patch. This
representation is the concatenation of the 2D coordinates of the $k$ vertices
together with the probability. The image coordinates are normalized, such that
their range is [-1,1].

---
Figure 3: The proposed method. Network $f$ contains an encoder-decoder network
with skip connections between them. Given an input Image I, it generates (1) a
low-resolution segmentation probability map $M$ (2) a patch-wise polygon array
$P$. Network $g$ contains (1) a neural renderer that employs a predefined
triangulation $T$ to generate the binary polygon segmentation map $R_{p}$ from
$P$ (2) an aggregation layer that rearranges $R_{p}$ according to the patch
location in the original image. The accumulated map is then gated, i.e.,
multiplied pixel-wise, by the initial probabilities provided by the upsampled
version of $M$ to generate the output segmentation map $M_{o}$.

A fixed triangulation $T$ is obtained by performing Delaunay triangulation to
a regular polygon of $k$ vertices. It is used to convert all polygons to a
raster image, assuming that locally the polygons are almost always close to
being convex. We set the z-axis to be constant and the perspective to be zero,
in order to use a 3D renderer for 2D objects.

The output $f(I)$ is divided into two parts, see Fig. 3. We denote by
$P=\\{P^{i}\\}$ the set of all $\frac{H}{s}\times\frac{W}{s}$ polygons, and by
$M$ the map of size $\frac{H}{s}\times\frac{W}{s}$ that contains all the
probabilities. The neural renderer $r$, given the vertices $P$ and the fixed
triangulation $T$, returns the polygon shape as a local image patch. In this
raster image, all pixels inside the polygon defined by $P$ through the
triangulation $T$ are assigned a value of one, zero for pixels outside, and
in-between values for pixels that the polygon divides:

${M^{i}_{p}}=r(P^{i},T)\in\mathbb{R}^{s\times s}\,,$ (1)

The aggregation layer collecting all local patches $M^{i}_{p}$, each in its
associated patch location, we obtain the map $R_{p}\in\mathbb{R}^{H\times W}$.
By performing nearest neighbor upsampling to the probability map $M$ by a
factor of $s$ in each dimension, and obtain a probability maps of the same
size as $R_{p}$. The final mask is the pixel-wise product between these two
maps.

${M}_{o}={\uparrow}M\odot{R_{p}}\,,$ (2)

where $\odot$ denotes elementwise multiplication and $\uparrow$ denotes
upsampling. This way, the low-resolution mask serves as a gating mechanism
that determines whether the polygon exists within an object of a certain class
in the output segmentation map.

To maximize the overlap between the output mask (before and after
incorporating the polygons) and the ground truth segmentation map $Y$, a soft
dice loss [31]:

$\mathcal{L}_{dice}(y,\hat{y})=1-\frac{2TP(y,\hat{y})+1}{2TP(y,\hat{y})+FN(y,\hat{y})+FP(y,\hat{y})+1}$
(3)

Where TP is the true positive between the ground truth $y$ and output mask
$\hat{y}$, FN is a false negative and FP is a false positive. In addition, a
binary cross-entropy is used for $M$:

$\mathcal{L}_{BCE}(y,\hat{y})=-y\log(\hat{y})+(1-y)\log(1-\hat{y})$ (4)

$\mathcal{L}=\mathcal{L}_{BCE}({\uparrow}M,Y)+\mathcal{L}_{dice}(M_{o},Y)$ (5)

This loss appears without weighting, in order to avoid an additional
parameter. During inference, we employ only $M_{o}$ and the second term can be
seen as auxiliary.

### Architecture

For the backbone $f_{1}$, the Harmonic Dense Net [4] is used, in which the
receptive field is of size $32\times 32$. The network contains six “HarD”
blocks, as described in [4], with 192, 256, 320, 480, 720, 1280 output
channels, respectively. The last block contains a dropout layer of 0.1. Pre-
trained ImageNet weights are used at initialization.

The decoder $f_{2}$ contains two upsampling blocks in order to obtain an
output resolution of $s=8$ times less than the original image size. Each block
contains two convolutional layers with a kernel size equal to 3 and zero
padding equal to one. In addition, we use batch normalization after the last
convolution layer before the activation function. The first layer’s activation
function is a ReLU, while the second layer’s activation functions are sigmoid
for the channel of the low-resolution map and $tanh$ for the polygon vertices.
Each layer receives a skip connection from the block of the encoder that has
the same spatial resolution. We note that our decoder requires considerably
fewer learnable parameters related to a regular decoder and, since only two
blocks are used, fewer skip connections from the encoder.

The subsequent network $g$ employs the mesh renderer of [16], with a zero
perspective, a camera distance of one, and output mask of size $s\times{s}$.
This way, the 3D renderer is adapted for the 2D rendering task.

In term of FLOPs, for image size of $256^{2}$ we get 16.51 GMacs for our model
and 25.44 GMacs for the comparable FCN. The peak memory consumption is 486MB
for our model and 506MB for the FCN.

Category | Full[train/val] | Occlusion[train/val]
---|---|---
Bicycle | 3501/1166 | 1029/385
Bus | 374/98 | 115/30
Car | 26226/4655 | 3333/914
Person | 16753/3395 | 1731/669
Train | 167/23 | 86/15
Truck | 481/93 | 108/19
Motorcycle | 707/149 | 175/55
Rider | 1711/544 | 291/155
Table 1: Cityscapes partitions training and validation set for each category, for the full dataset, and the occlusion subset. |  |  |
---|---|---|---
|  |  |
(a) | (b) | (c) | (d)
Figure 4: Sample segmentation masks on the Cityscapes validation dataset. (a) ground-truth (b) our method (c) FCN-HarDNet-85 [4] (d) the result of the active contour approach of [10]. The active contour method cannot handle occlusions or objects with overhanging parts. Method | Bike | Bus | Person | Train | Truck | Mcycle | Car | Rider | Mean
---|---|---|---|---|---|---|---|---|---
Polygon-RNN++ (with BS) [1] | 63.06 | 81.38 | 72.41 | 64.28 | 78.90 | 62.01 | 79.08 | 69.95 | 71.38
PSP-DeepLab [6] | 67.18 | 83.81 | 72.62 | 68.76 | 80.48 | 65.94 | 80.45 | 70.00 | 73.66
Polygon-GCN (with PS) [20] | 66.55 | 85.01 | 72.94 | 60.99 | 79.78 | 63.87 | 81.09 | 71.00 | 72.66
Spline-GCN (with PS) [20] | 67.36 | 85.43 | 73.72 | 64.40 | 80.22 | 64.86 | 81.88 | 71.73 | 73.70
Ours | 68.20 | 85.10 | 74.52 | 74.80 | 80.50 | 66.20 | 82.10 | 71.81 | 75.40
Deep contour [10] | 68.08 | 83.02 | 75.04 | 74.53 | 79.55 | 66.53 | 81.92 | 72.03 | 75.09
FCN-HarDNet-85 [4] | 68.26 | 84.98 | 74.51 | 76.60 | 80.20 | 66.25 | 82.36 | 71.57 | 75.59
Ours | 69.53 | 85.50 | 75.15 | 76.90 | 81.20 | 66.96 | 82.69 | 72.18 | 76.26
Table 2: Cityscapes segmentation results for two protocols: the top part refers to segmentation results with 15% expansion around the bounding box; the bottom part refers to segmentation results with no expansion around the bounding box Method | F1-Score | mIoU | WCov | FBound
---|---|---|---|---
FCN-UNet (Ronneberger 2015) | 87.40 | 78.60 | 81.80 | 40.20
FCN-ResNet34 | 91.76 | 87.20 | 88.55 | 75.12
FCN-HarDNet-85 (Chao 2019) | 93.97 | 88.95 | 93.60 | 80.20
DSAC (Marcos 2018) | - | 71.10 | 70.70 | 36.40
DarNet (Cheng 2019) | 93.66 | 88.20 | 88.10 | 75.90
Deep contour (Gur 2019) | 94.80 | 90.33 | 93.72 | 78.72
TDAC (Hatamizadeh 2020) | 94.26 | 89.16 | 90.54 | 78.12
Ours | 95.15 | 90.92 | 94.36 | 83.89
Table 3: Segmentation results on the Vaihingen building dataset. Figure 5: Performance under occlusion of our method (green diamond), in comparison to the leading UNet like network FCN-Hardnet85 [4] (blue circle) and the state-of-the-art active contour method [10] (red triangle). The x-axis denotes the amount of occlusion (see text). Shows is mean across all classes, see supplementary for a breakdown per class. |  |  |
---|---|---|---
|  |  |
(a) | (b) | (c) | (d)
Figure 6: Sample results of the proposed method on the building dataset. (a) The input image. (b) Ground truth segmentation. (c) The initial segmentation map $M$. (d) The final segmentation map $M_{o}$. First row: original images; second: a zoom-in view. Method | Dice | mIoU
---|---|---
FCN [2] | 28.84 | 28.71
U-Net [27] | 79.43 | 65.99
U-Net++[41] | 79.49 | 66.04
Res-UNet[37] | 79.49 | 66.07
Axial Attention U-Net[36] | 76.83 | 62.49
MedT[34] | 79.55 | 66.17
FCN-Hardnet85 | 79.52 | 66.06
Low res FCN-Hardnet85 ($M$) | 65.82 | 49.13
Ours ($M_{o}$) | 80.05 | 66.82
Table 4: Results for MoNu dataset |  |  |
---|---|---|---
|  |  |
|  |  |
(a) | (b) | (c) | (d)
Figure 7: Sample results of the proposed method on the MoNuSeg dataset. (a)
The input image. (b) Ground truth segmentation. (c) The low-resolution
segmentation map $M$. (d) The final segmentation map $M_{o}$.

##  4 Experiments

We present our results on three segmentation benchmarks and compare our
performance with the state-of-the-art methods. We further perform a study on
the effect of each parameter of the proposed algorithm, and on the effect of
the underlying network architecture. We put an added emphasis on the subset of
test images that depict an occlusion of the segmented object.

The Citycapes dataset [9] for instance segmentation includes eight object
classes with $5,000$ fine-annotated images, divided into $2,975$ images for
training, $500$ for validation, and $1,525$ for the test set.

The Cityscapes dataset contains urban scenes that are characterized by dense
and multiple objects. In order to fully understand the scene, the segmentation
method is required to be robust to occlusions of objects, which causes the
segmentation task to be more challenging. Following [20], we train and
evaluate given the bounding box of the object, such that the network outputs a
single mask for the entire instance.

To analyze the images with occlusions, the cases where the ground truth
contains at least two contours are identified. for this purpose, we define an
occlusion distance $d$ per image as:

$d=\max_{C_{i},C_{j}\in C}\min_{p\in C_{i},q\in C_{j}}\|p-q\|\,,$ (6)

where $C$ is the set of all contours of a specific class defined for a
specific image, each contour being a set of points. The distance $d$, which
degenerates to zero if there is only one contour, is then used to stratify the
test images by the number of occlusions. Tab. 1 contains the number of images
in both the validation set of the entire Cityscapes dataset, as well as for
the subset of which the occlusion distance $d$ is larger than one. As can be
seen, occlusions are fairly common for many of the objects.

The baseline methods we compare to on the Cityscapes dataset are
PSPDeepLab[6], Polygon-RNN++[1], Curve-GCN[20] and Deep active contours [10].
Another baseline, denoted as “FCN-HarDNet-85” employs a fully-convolutional
network with the same backbone [4] we used for our own method.

We note that our method employs 17M parameters in its decoder, while the “FCN-
HarDNet-85” employs a decoder with 22M parameters.

The Vaihingen [28] dataset consists of 168 aerial images, at a resolution of
$256\times 256$, from this German city. The task is to segment the central
building in a very dense and challenging environment, which includes other
structures, streets, trees, and cars. The dataset is divided into 100
buildings for training, and the remaining 68 for testing.

For building segmentation, we use an encoder-decoder architecture based on the
HarDNet-85 [4] backbone and produce $M$ at 1/8 resolution ($s=8$).

The baseline methods we compare to on the Vaihingen buildings dataset are
“DSAC”[23], “DarNet”[7], “TDAC”[11] and Deep active contours [10]. We also
present results for “FCN-UNET” [27], “FCN-ResNet-34” and “FCN-HarDNet-85”, all
based on fully-convolutional network decoders.

The MoNuSeg dataset [18] contains a training set with 30 microscopic images
from seven organs with annotations of 21,623 individual nuclei, the test
dataset contains 14 similar images. Following previous work, we resized the
images into a resolution of $512\times 512$ [34]. An encoder-decoder
architecture based on the HarDNet-85 [4] backbone. Since the recovered
elements are very small, our method employs a map $M$ at 1/4 resolution
($s=4$), and not lower resolutions as in other datasets.

The baseline methods we compare to on the Monu dataset are “FCN”[2],
“UNET”[27], “UNET++”[41], Res-Unet[37], Axial attention Unet[36] and Medical
transformer [34].

Finally, we report results on CUB-200-2011 [35]. The dataset contains 5994
images of different birds for training, and 5794 for validation. We trained
the model for the same setting of cityscapes where the input image size is
$256^{2}$ and the evaluation is on the original image size.

Training details  For training our network, we use the ADAM optimizer with an
initial learning rate of 0.0003 and a gamma decay of 0.7, applied every 50
epochs. The batch size is 32 and the parameter of the weight decay
regularization is set to $5\cdot 10^{-5}$. A single GeForce RTX 2080 Ti is
used for training on all datasets. Our networks are trained for up to 1000
epochs.

Following [10], we used for the Cityscapes images a set of augmentations that
includes: (i) scaling by a random scale factor in the range of $[0.75,1.25]$
and then cropping an image of size $256\times 256$, (ii) color jitter with the
parameters of brightness sampled uniformly between $[0,0.4]$, a contrast in
the range $[0,0.4]$, saturation in the range $[0,0.4]$, and hue in the range
$[0,0.1]$, (iii) a random horizontal flip with a probability of 0.5.

For the Vaihingen dataset, also following [10], we used (i) a random rotation
augmentation of 360 degrees and a scale range of $[0.75,1.5]$, (ii) a random
vertical and horizontal flip with a probability of 0.5, (iii) random color
jitter with a maximal value of 0.6 for brightness, 0.5 for contrast, 0.4 for
saturation, and 0.025 for hue.

For the MoNu dataset, we used (i) a random rotation augmentation of ±20
degrees and a scale range of $[0.75,1.25]$, (ii) a random horizontal flip with
a probability of 0.5, (iii) random color jitter with a maximal value of 0.4
for brightness, 0.4 for contrast, 0.4 for saturation, and 0.1 for hue.

Evaluation Metrics  The Cityscapes segmentation results employ the common
metrics of mean Intersection-over-Union (IoU). The Vaihingen benchmark is
equipped with a set of additional evaluation metrics. The results are
evaluation by the IoU of the single class, and the F1-score of retrieving
building pixels. The weighted Coverage (WCov) score is computed as the IoU for
each object in the image weighted by the related size of the object from the
whole image. Finally, the Boundary F-score (BoundF) is the average of the
F1-scores computed with thresholds ranging from 1 to 5 pixels around the
ground truth boundaries, as described by [28].

Results  The Cityscapes results are presented in Tab. 2. As can be seen, our
method outperforms all baseline methods, including active contour methods and
encoder-decoder architectures. Furthermore, one can observe the gap in
performance between our method and the U-Net-like method that is based on the
same backbone architecture (“FCN-HarDNet-85”). This is despite the reduction
in the number of trainable parameters and the size of the representation (each
$8\times 8$ patch is represented by $k+1$ output floats, instead of 64
floats).

Fig. 4 presents sample results and compares our method to the leading active
contour method (“Deep contour”) and to the FCN-HarDNet-85 baseline. As can be
seen, the active contour method struggles to model objects that deviate
considerably from their convex hull, and cannot handle occlusions well. In
comparison to the FCN baseline, our method provides more accurate details
(consider the bicycle fork in the first row), but the overall shape is the
same.

The results for images stratified by the level of occlusion (the measurement
$d$) are presented in Fig. 5. All methods present a decrease in performance as
the amount of occlusion increases. The FCN-HarDNet-85 maintains a stable
performance level below our method, and the active contour method suffers the
most from occlusion, as expected.

The results for the Vaihingen building dataset are presented in Tab. 3. Our
method outperforms both the fully connected segmentation methods, as well as
the active-contour based methods. Fig. 6 presents typical samples,
demonstrating the refinement that is obtained by the polygons over the initial
segmentation mask.

The results for MoNu dataset are reported in Tab. 4. We outperform all
baselines for both the Dice score and mean-IoU. Our algorithm also obtains
better performances in comparison to the fully convolutional segmentation
network with the same backbone Hardnet-85. The improvement from low resolution
mask ($M$) to the output mask ($M_{o}$) is from 49.13% IoU to 66.82%. The full
resolution FCN-Hardnet85 obtains a lower result of 66.06 IoU. Fig. 7 presents
the output mask of our algorithm for samples from the test set.

The results for CUB-200-2011 are reported in Tab. 5. We obtain better results
than both FCN-Hardnet85 and Deep Contour on all metrics. Fig. 8 present the
output mask of our algorithm for samples from CUB-200-2011.

Method | Dice | mIoU | FBound
---|---|---|---
U-Net [27] | 93.84 | 88.80 | 76.54
Deep contour [10] | 93.72 | 88.35 | 75.92
FCN-Hardnet85 | 94.91 | 90.51 | 79.73
Low res FCN-Hardnet85 ($M$) | 84.45 | 73.43 | 39.47
Ours ($M_{o}$) | 95.11 | 90.91 | 79.86
Table 5: Results for CUB dataset |  |  |
---|---|---|---
Figure 8: Sample results of our method on the CUB dataset Figure 9: The effect
on accuracy of varying the number of vertices to be one of [4, 5, 8, 16, 24]
and the patch size to be in the set [4, 8, 16, 32]. Results are reported on
the Vaihingen building dataset.

Parameter sensitivity  The proposed patch-wise polygon scheme adds two
hyperparameters: the size of the patch that each polygon for each class is
applied to and the number of vertices $k$ in each polygon. In order to study
the effect of these parameters, we varied the number of vertices to be in the
range of $[4,5,8,16,24]$ and the size of the patch to be one of $[4\times
4,8\times 8,16\times 16,32\times 32]$.

Fig. 9 shows our network performances for building segmentation, when varying
the two parameters. As can be seen, for small patches, for instance, $4\times
4$ or $8\times 8$, the number of nodes that optimizes performance is five. As
may be expected, for larger patches, the method performs better with a higher
number of vertices.

We also checked our performances for different Hardnet backbones -
HarDNnet39DS, HarDNnet68DS, HarDNnet68, HarDNnet85 (DS stands for for depth-
wise-separable, see [4]). The results appear in Tab. 6 for the building
dataset. As expected, the performance improves with the increase in the number
of parameters in the network. It can also be observed that the low- resolution
mask after upsampling $M$ and the refined mask $M_{o}$ have a performance gap
of 4% in terms of mean IoU, in favor of the latter, across all of the
backbones.

Backbone | Params[#] |  $M_{o}$ [mIoU] |  $M$ [mIoU]
---|---|---|---
HarDNet-39DS | 3.5M | 89.11 | 85.07
HarDNet-68DS | 4.2M | 89.75 | 85.40
HarDNet-68 | 17.6M | 90.50 | 85.60
HarDNet-85 | 36.7M | 90.92 | 85.98
Table 6: Backbone comparison for the building segmentation task using various
HarDNet variants [4]. Showing the results of the upscaled gating signal $M$
and the final output $M_{0}$ that includes the polygons’ boundaries.

##  5 Discussion

The ability to work with polygons allows our network to produce segmentation
maps at a resolution that is limited only by machine precision. This can be
seen in Fig. 2(b), where the polygons have fractional coordinates. These
polygons can be rasterized at any resolution.

The property of practically-infinite resolution may be of use in applications
that require an adaptive resolution, depending on the object or its specific
regions. An example is foreground segmentation for background replacement
(virtual “green screens”), in which hair regions require finer resolutions.
Such scenarios would raise research questions about training networks at a
resolution that is higher than the resolution of the ground truth data, and on
the proper evaluation of the obtained results.

We note that while we do not enforce the polygons of nearby patches to be
compatible at the edges, this happens naturally, as can be seen in Fig. 2(b).
Earlier during development, we have attempted to use loss terms that encourage
such compatibility but did not observe an improvement.

Further inspection of the obtained polygons reveals that polygons with $k>4$
often employ overlapping vertices in order to match the square cell edges. In
empty regions, polygons tend to become of zero-area by sticking to one of the
boundaries. These behaviors emerge without being directly enforced. In some
cases, phantom objects are detected by single polygons, as in the top right
part of Fig. 2(b). These polygons are removed by the gating process that
multiplies with $M$, as can be seen in Fig. 2(c). We do not observe such
polygons in the output images.

##  6 Conclusions

The pixel grid representation is commonly used by deep segmentation networks.
In this work, we present an alternative approach that encodes the local
segmentation maps as polygons. The method employs a neural renderer to
backprop the rasterization process. A direct comparison to a method that
employs the same backbone but without the polygonal representation reveals a
significant improvement in performance on multiple benchmarks, despite using
fewer training parameters. Moreover, our method also outperforms all other
segmentation methods.

## Acknowledgment

This project has received funding from the European Research Council (ERC)
under the European Unions Horizon 2020 research and innovation programme
(grant ERC CoG 725974).

## References

  * [1] David Acuna, Huan Ling, Amlan Kar, and Sanja Fidler.  Efficient interactive annotation of segmentation datasets with polygon-rnn++.  In Proceedings of the IEEE conference on Computer Vision and Pattern Recognition, pages 859–868, 2018.
  * [2] Vijay Badrinarayanan, Alex Kendall, and Roberto Cipolla.  Segnet: A deep convolutional encoder-decoder architecture for image segmentation.  IEEE transactions on pattern analysis and machine intelligence, 39(12):2481–2495, 2017.
  * [3] Vicent Caselles, Ron Kimmel, and Guillermo Sapiro.  Geodesic active contours.  International journal of computer vision, 22(1):61–79, 1997.
  * [4] Ping Chao, Chao-Yang Kao, Yu-Shan Ruan, Chien-Hsiang Huang, and Youn-Long Lin.  Hardnet: A low memory traffic network.  In Proceedings of the IEEE International Conference on Computer Vision, pages 3552–3561, 2019.
  * [5] Liang-Chieh Chen, George Papandreou, Iasonas Kokkinos, Kevin Murphy, and Alan L Yuille.  Semantic image segmentation with deep convolutional nets and fully connected crfs.  arXiv preprint arXiv:1412.7062, 2014.
  * [6] Liang-Chieh Chen, George Papandreou, Iasonas Kokkinos, Kevin Murphy, and Alan L Yuille.  Deeplab: Semantic image segmentation with deep convolutional nets, atrous convolution, and fully connected crfs.  IEEE transactions on pattern analysis and machine intelligence, 40(4):834–848, 2017.
  * [7] Dominic Cheng, Renjie Liao, Sanja Fidler, and Raquel Urtasun.  Darnet: Deep active ray network for building segmentation.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 7431–7439, 2019.
  * [8] Laurent D Cohen.  On active contour models and balloons.  CVGIP: Image understanding, 53(2):211–218, 1991.
  * [9] Marius Cordts, Mohamed Omran, Sebastian Ramos, Timo Rehfeld, Markus Enzweiler, Rodrigo Benenson, Uwe Franke, Stefan Roth, and Bernt Schiele.  The cityscapes dataset for semantic urban scene understanding.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 3213–3223, 2016.
  * [10] Shir Gur, Tal Shaharabany, and Lior Wolf.  End to end trainable active contours via differentiable rendering.  arXiv preprint arXiv:1912.00367, 2019.
  * [11] Ali Hatamizadeh, Debleena Sengupta, and Demetri Terzopoulos.  End-to-end trainable deep active contour models for automated image segmentation: Delineating buildings in aerial imagery.  In European Conference on Computer Vision, pages 730–746. Springer, 2020.
  * [12] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [13] Namdar Homayounfar, Wei-Chiu Ma, Shrinidhi Kowshika Lakshmikanth, and Raquel Urtasun.  Hierarchical recurrent attention networks for structured online maps.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 3417–3426, 2018.
  * [14] Gao Huang, Zhuang Liu, Laurens Van Der Maaten, and Kilian Q Weinberger.  Densely connected convolutional networks.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 4700–4708, 2017.
  * [15] Michael Kass, Andrew Witkin, and Demetri Terzopoulos.  Snakes: Active contour models.  International journal of computer vision, 1(4):321–331, 1988.
  * [16] Hiroharu Kato, Yoshitaka Ushiku, and Tatsuya Harada.  Neural 3d mesh renderer.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 3907–3916, 2018.
  * [17] Satyanad Kichenassamy, Arun Kumar, Peter Olver, Allen Tannenbaum, and Anthony Yezzi.  Gradient flows and geometric active contour models.  In Proceedings of IEEE International Conference on Computer Vision, pages 810–815. IEEE, 1995.
  * [18] Neeraj Kumar, Ruchika Verma, Deepak Anand, Yanning Zhou, Omer Fahri Onder, Efstratios Tsougenis, Hao Chen, Pheng-Ann Heng, Jiahui Li, Zhiqiang Hu, et al.  A multi-organ nucleus segmentation challenge.  IEEE transactions on medical imaging, 39(5):1380–1391, 2019.
  * [19] Justin Liang, Namdar Homayounfar, Wei-Chiu Ma, Yuwen Xiong, Rui Hu, and Raquel Urtasun.  Polytransform: Deep polygon transformer for instance segmentation.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 9131–9140, 2020.
  * [20] Huan Ling, Jun Gao, Amlan Kar, Wenzheng Chen, and Sanja Fidler.  Fast interactive object annotation with curve-gcn.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 5257–5266, 2019.
  * [21] Jonathan Long, Evan Shelhamer, and Trevor Darrell.  Fully convolutional networks for semantic segmentation.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 3431–3440, 2015.
  * [22] Matthew M Loper and Michael J Black.  Opendr: An approximate differentiable renderer.  In European Conference on Computer Vision, pages 154–169. Springer, 2014.
  * [23] Diego Marcos, Devis Tuia, Benjamin Kellenberger, Lisa Zhang, Min Bai, Renjie Liao, and Raquel Urtasun.  Learning deep structured active contours end-to-end.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 8877–8885, 2018.
  * [24] Gabriel L Oliveira, Claas Bollen, Wolfram Burgard, and Thomas Brox.  Efficient and robust deep networks for semantic segmentation.  The International Journal of Robotics Research, 37(4-5):472–491, 2018.
  * [25] Vivek Parmar, Narayani Bhatia, Shubham Negi, and Manan Suri.  Exploration of optimized semantic segmentation architectures for edge-deployment on drones.  arXiv preprint arXiv:2007.02839, 2020.
  * [26] Chao Peng, Xiangyu Zhang, Gang Yu, Guiming Luo, and Jian Sun.  Large kernel matters–improve semantic segmentation by global convolutional network.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 4353–4361, 2017.
  * [27] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In International Conference on Medical image computing and computer-assisted intervention, pages 234–241. Springer, 2015.
  * [28] Franz Rottensteiner, Gunho Sohn, Jaewook Jung, Markus Gerke, Caroline Baillard, Sébastien Bénitez, and U Breitkopf.  International society for photogrammetry and remote sensing, 2d semantic labeling contest.  http://www2.isprs.org/commissions/comm3/wg4/semantic-labeling.html.
  * [29] Christian Rupprecht, Elizabeth Huaroc, Maximilian Baust, and Nassir Navab.  Deep active contours.  arXiv preprint arXiv:1607.05074, 2016.
  * [30] Karen Simonyan and Andrew Zisserman.  Very deep convolutional networks for large-scale image recognition.  arXiv preprint arXiv:1409.1556, 2014.
  * [31] Carole H Sudre, Wenqi Li, Tom Vercauteren, Sebastien Ourselin, and M Jorge Cardoso.  Generalised dice overlap as a deep learning loss function for highly unbalanced segmentations.  In Deep learning in medical image analysis and multimodal learning for clinical decision support, pages 240–248. Springer, 2017.
  * [32] Christian Szegedy, Wei Liu, Yangqing Jia, Pierre Sermanet, Scott Reed, Dragomir Anguelov, Dumitru Erhan, Vincent Vanhoucke, and Andrew Rabinovich.  Going deeper with convolutions.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 1–9, 2015.
  * [33] Michael Treml, José Arjona-Medina, Thomas Unterthiner, Rupesh Durgesh, Felix Friedmann, Peter Schuberth, Andreas Mayr, Martin Heusel, Markus Hofmarcher, Michael Widrich, et al.  Speeding up semantic segmentation for autonomous driving.  In MLITS, NIPS Workshop, volume 2, 2016.
  * [34] Jeya Maria Jose Valanarasu, Poojan Oza, Ilker Hacihaliloglu, and Vishal M Patel.  Medical transformer: Gated axial-attention for medical image segmentation.  arXiv preprint arXiv:2102.10662, 2021.
  * [35] Catherine Wah, Steve Branson, Peter Welinder, Pietro Perona, and Serge Belongie.  The Caltech-UCSD birds-200-2011 dataset.  2011\.
  * [36] Huiyu Wang, Yukun Zhu, Bradley Green, Hartwig Adam, Alan Yuille, and Liang-Chieh Chen.  Axial-deeplab: Stand-alone axial-attention for panoptic segmentation.  In European Conference on Computer Vision, pages 108–126. Springer, 2020.
  * [37] Xiao Xiao, Shen Lian, Zhiming Luo, and Shaozi Li.  Weighted res-unet for high-quality retina vessel segmentation.  In 2018 9th international conference on information technology in medicine and education (ITME), pages 327–331. IEEE, 2018.
  * [38] Fisher Yu and Vladlen Koltun.  Multi-scale context aggregation by dilated convolutions.  arXiv preprint arXiv:1511.07122, 2015.
  * [39] Hang Zhang, Kristin Dana, Jianping Shi, Zhongyue Zhang, Xiaogang Wang, Ambrish Tyagi, and Amit Agrawal.  Context encoding for semantic segmentation.  In Proceedings of the IEEE conference on Computer Vision and Pattern Recognition, pages 7151–7160, 2018.
  * [40] Hengshuang Zhao, Jianping Shi, Xiaojuan Qi, Xiaogang Wang, and Jiaya Jia.  Pyramid scene parsing network.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2881–2890, 2017.
  * [41] Zongwei Zhou, Md Mahfuzur Rahman Siddiquee, Nima Tajbakhsh, and Jianming Liang.  Unet++: A nested u-net architecture for medical image segmentation.  In Deep learning in medical image analysis and multimodal learning for clinical decision support, pages 3–11. Springer, 2018.


#  16O spectral function from coupled-cluster theory:
applications to lepton-nucleus scattering

J. E. Sobczyk  Institut für Kernphysik and PRISMA+ Cluster of Excellence,
Johannes Gutenberg-Universität, 55128 Mainz, Germany  S. Bacca  Institut für
Kernphysik and PRISMA+ Cluster of Excellence, Johannes Gutenberg-Universität,
55128 Mainz, Germany  Helmholtz-Institut Mainz, Johannes Gutenberg-Universität
Mainz, D-55099 Mainz, Germany

###### Abstract

We calculate the 16O spectral function by combining coupled-cluster theory
with a Gaussian integral transform and by expanding the integral kernel in
terms of Chebyshev polynomials to allow for a quantification of the
theoretical uncertainties. We perform an analysis of the spectral function and
employ it to predict lepton-nucleus scattering. Our results well describe the
16O electron scattering data in the quasi-elastic peak for momentum transfers
$|\mathbf{q}|\gtrapprox 500$ MeV and electron energies up to 1.2 GeV,
extending therefore the so-called first principles approach to lepton-nucleus
cross sections well into the relativistic regime. To prove the applicability
of this method to neutrino-nucleus cross sections, we implement our 16O
spectral functions in the NuWro Monte Carlo event generator and provide a
comparison with recently published T2K neutrino data.

##  I Introduction

Long-baseline neutrino-oscillation programs hold a special place among
neutrino experimental endeavours. Driven by their ambitious aims—measuring the
CP-violating phase in the electroweak sector and searching for new physics—the
next-generation experiments DUNE Acciarri  _et al._ (2015) and T2HK Abe  _et
al._ (2015) will drastically reduce the statistical uncertainty to keep them
at the level of a few percent. This implies that systematical uncertainties,
until now overshadowed by the statistics, will require more attention. A
considerable source of uncertainty comes from the modeling of neutrino-nucleus
interactions Alvarez-Ruso  _et al._ (2018); Ankowski  _et al._ (2022); Ruso
_et al._ (2022). Given that neutrinos are elusive particles, they can in fact
be detected only indirectly via the measurement of the final particles
produced in their interaction with target nuclei. Moreover, oscillation
experiments use neutrino fluxes of a wide energy-range, from hundreds of MeV
to few GeV, which make them sensitive to a variety of dynamical mechanisms.
These conditions require an excellent understanding of the underlying
processes to precisely reconstruct the neutrino energy in each observed event.

A theoretical description of electroweak reactions that starts from the forces
among nucleons and their interactions with external probes, and is based on a
numerical solution of the problem within controlled approximations, is
arguably the doorway to a solid understanding of the dynamical mechanism
governing lepton-nucleus scattering. This so-called first principles (or ab
initio) approach is on the one hand computationally intensive, but on the
other hand offers the prospects of quantifying and possibly reducing nuclear
physics uncertainties in the computed cross sections. For light nuclei, Green
Function Monte Carlo (GFMC) has been very successful in delivering predictions
of electroweak cross-sections Lovato  _et al._ (2018, 2020). For nuclei up to
mass number 40, the Lorentz integral transform combined with coupled-cluster
theory (LIT-CC) has been recently proven to work well Sobczyk  _et al._ (2020,
2021). Both the above mentioned methods are capable of describing the low
energy-momentum part of the lepton-nucleus cross section which depends on the
details of the nuclear dynamics in the final state, the so-called final state
interaction (FSI). In particular, in the case of the LIT-CC method an
extension to higher energies and momenta is complicated by the necessity of
using soft nuclear Hamiltonians which typically have a cutoff of about 500 MeV
and therefore cannot be reliably used to describe FSI beyond that point.
However, at higher energies while ground-state correlations remain important,
FSIs become negligible and one can assume that only one nucleon in the nucleus
interacts with the external probe and is knocked out after getting all the
momentum transferred from the lepton. In this regime, one can use the spectral
function (SF) to compute cross sections, leading to a simplification of the
computational task.

The SF formalism is a well established approach to describe lepton-nucleus
scattering within the impulse approximation (IA). It is based on a
factorization ansatz of the ground-state nucleus in terms of one nucleon which
participates to the interaction vertex with the external probe, while the
remaining $A-1$ nucleus is a spectator. The SF formalism is amenable to an
extension of the first principle description to the relativistic regime,
because one can use relativistic currents in the interaction vertex and even
account for higher-energy mechanisms like the pion production Rocco  _et al._
(2019). This is particularly useful in neutrino physics, allowing to address
consistently various reactions within the same underlying formalism.

The experimental collaborations perform neutrino energy reconstruction using
Monte Carlo (MC) event generators. The analysis requires the knowledge of the
semi-inclusive reactions, most importantly the outgoing protons and pions
distributions. Currently, the neutrino community is devoting considerable
efforts into improving the quality of implemented nuclear models and into
going beyond the simple inclusive cross-sections. Also in this respect the SF
formalism is a convenient tool to model the electroweak processes, because it
allows to address semi-inclusive knockout reactions in a straightforward way.
For example, phenomenological SFs can be constructed using the experimental
$(e,e^{\prime}p)$ data, see, e.g., Ref. Jiang  _et al._ (2022) for the recent
results from the Jefferson Laboratory. There are several theoretical models of
phenomenological nature of the spectral function already available for
neutrino studies Benhar  _et al._ (1994); Nieves and Sobczyk (2017); Buss  _et
al._ (2007). However, with the prospect of quantifying theoretical
uncertainties, it is worth investing into the development of spectral
functions derived from first principles. A first calculation based on the self
consistent Green’s function (SCGF) method was provided in Refs. Rocco and
Barbieri (2018); Barbieri  _et al._ (2019). Recently, a new method to
construct spectral functions for the many-body system based on Chebyshev
polynomials expansion of the integral kernel (ChEK method) Roggero (2020);
Sobczyk and Roggero (2021) was proposed. The ChEK method used in conjunction
with coupled-cluster theory Hagen  _et al._ (2014) was benchmarked on 4He
leading to a good agreement with electron-scattering data Sobczyk  _et al._
(2022). In this paper, we present a computation of the 16O spectral function
with the ChEK method. The main advantage of this approach is that it accounts
for the uncertainties of the spectral reconstruction. This feature becomes
especially valuable when theoretical uncertainties are propagated to the
computed cross sections, allowing for a comparison to experimental data on
equal footing.

The paper is organized in the following way. In Sec. II we review the quasi-
elastic (QE) process within impulse approximation. In Sec. III we present the
theoretical framework in which we perform the calculation of the spectral
functions. In Sec. IV we apply the formalism to 16O, comparing our results
both to electron and neutrino scattering data. The latter is done withing the
NuWro Monte Carlo generator Juszczak  _et al._ (2006); Golan  _et al._ (2012).
Finally, we conclude in Sec. V.

##  II Spectral function formalism in quasi-elastic scattering

The differential cross-section for the lepton-nucleus scattering can be
expressed as

$\frac{d^{2}\sigma}{d\omega
d\cos\theta}=\kappa\frac{|\mathbf{k}|}{|\mathbf{k}^{\prime}|}L_{\mu\nu}W^{\mu\nu}\,,$
(1)

where the energy-momentum transfer is given by $q=(\omega,\mathbf{q})$, the
scattering angle $\theta$ and initial and final lepton four-momenta are
$k=(E_{k},\mathbf{k})$ and $k^{\prime}=(E_{k^{\prime}},\mathbf{k^{\prime}})$,
respectively. The interaction vertex depends on the process with

$\kappa_{\mathrm{EM}}=\left(\frac{\alpha}{q^{2}}\right)^{2}\,,~{}~{}\kappa_{\mathrm{CC,NC}}=\left(\frac{G_{F}\cos\theta_{C}}{2\pi}\right)^{2}\,,$
(2)

for electromagnetic (EM), charge-current (CC) or neutral-current (NC),
respectively. The lepton tensor is given by

$L_{\mu\nu}=2a[k_{\mu}k_{\nu}^{\prime}+k_{\mu}^{\prime}k_{\nu}-g_{\mu\nu}(kk^{\prime})\pm
i\eta\epsilon_{\mu\nu\alpha\beta}k^{{}^{\prime}\alpha}k^{\beta}]\,,$ (3)

with $a=1$, $\eta=0$ for electromagnetic and $a=4$, $\eta=1$ for electroweak
reactions. The hadronic tensor

$W^{\mu\nu}=\sum_{f}\delta^{4}(p_{0}+q-p_{f})\langle
0|\left(J^{\mu}\right)^{\dagger}|\Phi_{f}\rangle\langle\Phi_{f}|J^{\nu}|0\rangle\,,$
(4)

where $J_{\mu}$ is the electroweak current, and $|0\rangle$ and
$|\Phi_{f}\rangle$ are the initial and final nuclear state with respective
four momenta $p_{0}$ and $p_{f}$, depends on the reaction mechanism under
examination. For the electroweak current $J_{\mu}$, we consider here only one-
body operators which in the notation of the second quantization can be written
as

$J^{\mu}=\sum_{\alpha,\beta}\langle\beta|j^{\mu}|\alpha\rangle
a_{\beta}^{\dagger}a_{\alpha}\,,$ (5)

where $\alpha$ and $\beta$ are the quantum numbers of single-particle states
created and annihilated by the respective operators $a^{\dagger}$ and $a$.
Within the spectral function formalism, we can use the fully relativistic
currents in the matrix element

$\langle p+q|j^{\mu}|p\rangle=\bar{u}(p+q)\left(V^{\mu}+A^{\mu}\right)u(p)\,,$
(6)

with Dirac spinors $u$ and the single-nucleon current $j^{\mu}$ having a
vector-axial structure. Constructing the most general form of $V^{\mu}$ and
$A^{\mu}$ using the available four-vectors, we have

$\begin{split}V^{\mu}=F_{1}\gamma^{\mu}+\frac{F_{2}}{2m}i\sigma^{\mu\nu}q_{\nu}\,,\\\
A^{\mu}=F_{A}\gamma^{\mu}\gamma^{5}+\frac{F_{P}}{m}q^{\mu}\gamma^{5}\,,\end{split}$
(7)

where form factors denoted by $F$ depend on the considered process. For the EM
scattering we will use $F_{1}^{n,p}$, $F_{2}^{n,p}$ parametrized as in Ref.
Bradford  _et al._ (2006). The CC vector form-factor is related to the
electromagnetic ones as $F_{i}=F_{i}^{p}-F_{i}^{n}$. The axial form factors –
present only in the weak interactions – are related under PCAC (partially
conserved axial current):

$F_{P}(Q^{2})=\frac{2m^{2}}{Q^{2}+m_{\pi}^{2}}{F_{A}(Q^{2})}\,,$ (8)

with $F_{A}$ taken as a dipole with $M_{A}=1030$ MeV axial mass.

Under the assumption that the struck nucleon does not interact with the final
nuclear state, we can factorize the final plane-wave nucleon with momentum
$\mathbf{p}^{\prime}$ and the $A-1$ nuclear state as $|\Phi_{f}\rangle\to
a_{p^{\prime}}|\Phi_{A-1}\rangle$. By inserting a complete set of intermediate
states, $\int d^{3}\mathbf{p}/(2\pi)^{3}|p\rangle\langle
p|a_{p}a^{\dagger}_{p}$, the many-body matrix element of the current operator
can be approximated by

$\begin{split}\langle\Phi_{f}|J^{\nu}|0\rangle\approx\int\frac{d^{3}\mathbf{p}}{(2\pi)^{3}}\langle
p^{\prime}|j^{\mu}|p\rangle\sum_{\alpha}\langle\mathbf{p}|\alpha\rangle\langle\Phi_{A-1}|a_{\alpha}|0\rangle\,.\end{split}$
(9)

The hadron tensor factorizes the interaction vertex while the nuclear effects
are encapsulated into the spectral function $S(\mathbf{p},E)$ 111To simplify
the notation, from now on we will write $\mathbf{p}$ instead of $|\mathbf{p}|$
when referring to an argument of the SF, $S(\mathbf{p},E)$, and momentum
distribution, $n(\mathbf{p})$. , separately for neutrons and protons, as

$\displaystyle W^{\mu\nu}(q)$
$\displaystyle=\int\frac{d^{3}\mathbf{p}}{(2\pi)^{3}}dE\frac{m}{E_{p}}\frac{m}{E_{p+q}}$
$\displaystyle\big{[}S^{n}(\mathbf{p},E)w_{n}^{\mu\nu}(p,q)+S^{p}(\mathbf{p},E)w_{p}^{\mu\nu}(p,q)\big{]}$
$\displaystyle\times\delta(\omega+E-E_{p+q}-E_{f}^{kin})\,,$ (10)

where $E_{f}^{kin}$ is the kinetic energy of $A-1$ system and $E_{p+q}$ is the
kinetic energy of outgoing nucleon. Here, the spectral function is defined as

$\begin{split}&S(\mathbf{p},E)=\sum_{\alpha,\beta}\langle\mathbf{p}|\alpha\rangle\langle\mathbf{p}|\beta\rangle^{\dagger}\\\
&\sumint_{\Phi_{A-1}}\langle
0|a_{\beta}^{\dagger}|\Phi_{A-1}\rangle\langle\Phi_{A-1}|a_{\alpha}|0\rangle\delta\big{(}E-(E_{0}-E_{\Phi})\big{)}\,.\end{split}$
(11)

and it gives the probability distribution of kicking a nucleon with momentum
$|\mathbf{p}|$ out of the ground-state, leaving it with an excitation energy
$E$. Finally, the factorized interaction vertex $w^{\mu\nu}(p,q)=\langle
p+q|j^{\mu}|p\rangle^{\dagger}\langle p+q|j^{\nu}|p\rangle$ employs the
current from Eq. (6).

Within the IA, the outgoing nucleon is decoupled from the nuclear final state,
hence one is able to factorize the high-energy physics taking place at the
interaction vertex from the properties of the nuclear ground state. However,
it is well known that neglecting the FSI at the intermediate momentum transfer
of the order of hundreds of MeV leads to some inconsistencies with the data
Sobczyk (2017). The QE peak is shifted towards higher energy transfers and it
exhibits too much strength. These inconsistencies can be partially alleviated
if the FSI for the struck nucleon are also included by introducing an optical
potential to describe the interaction of the struck nucleon with the rest of
the nucleus. Its real part amounts to the potential energy of the nucleon in
the nuclear medium, while the imaginary part is responsible to account for the
absorption channels of outgoing nucleon.

Phenomenological optical potentials are fitted using elastic nucleon-nucleus
scattering data. In the present calculations we will use the real part of the
optical potential for 16O from Ref. Cooper  _et al._ (1993). The relativistic
potential gives the scalar and vector contributions to the Dirac equation,
dependent on the kinetic energy and the radial position. Following the same
steps as in Ref. Ankowski  _et al._ (2015), averaging over the density of
protons, we arrive at the real part of optical potential
$\mathrm{Re}\mathrm{U}(t_{\mathrm{kin}})$ shown in Fig. 1 with
$t_{\mathrm{kin}}$ being the kinetic energy of the outgoing nucleon. The fit
of Ref. Cooper  _et al._ (1993) is reliable only above $t_{\mathrm{kin}}=25$
MeV (momentum $|\mathbf{p}^{\prime}|\approx 215$ MeV). The inclusion of the
imaginary part introduced as a folding of the cross section with a Lorentzian
function Ankowski  _et al._ (2015), was reported to overestimate the
absorption rate leading to non-physical large tails from the Lorentzian
distribution, see Ref. Benhar  _et al._ (1991). Therefore, we will not take
the imaginary part into account in our current predictions, since the topic
requires further investigations.

In essence, within our treatment of the FSI, the real part of optical
potential enters the hadron tensor changing the energy conservation as

$\displaystyle W^{\mu\nu}_{\mathrm{FSI}}(q)$
$\displaystyle=\int\frac{d^{3}\mathbf{p}}{(2\pi)^{3}}dE\frac{m}{E_{p}}\frac{m}{E_{p+q}}$
$\displaystyle\big{[}S^{n}(\mathbf{p},E)w_{n}^{\mu\nu}(p,q)+S^{p}(\mathbf{p},E)w_{p}^{\mu\nu}(p,q)\big{]}$
$\displaystyle\times\delta(\omega+E-E_{p+q}-E_{f}^{kin}-\mathrm{Re}\mathrm{U})\,.$
(12)

Following Ref. Ankowski  _et al._ (2015), we take $\mathrm{Re}\mathrm{U}$ at
$t_{kin}=\sqrt{m^{2}+\mathbf{q}^{2}}-m^{2}$.

Figure 1: The real part of optical potential for 16O as a function of the
kinetic energy of the outgoing proton from Ref. Cooper  _et al._ (1993).

##  III Spectral functions from the ChEK method

Spectral functions as shown in Eq. (11) are defined in terms of the imaginary
part of the hole propagator in a many-body system

$\begin{split}&\mathrm{Im}G_{h}(\alpha,\beta,E)=\\\
&-\pi\sumint_{\Phi_{A-1}}\langle
0|a_{\beta}^{\dagger}|\Phi_{A-1}\rangle\langle\Phi_{A-1}|a_{\alpha}|0\rangle\delta\big{(}E-(E_{0}-E_{\Phi})\big{)}\,.\end{split}$
(13)

The spectrum of excited states $\Phi_{A-1}$ contains bound and continuum
states, making the direct calculation of $\mathrm{Im}G_{h}(\alpha,\beta,E)$
challenging. To circumvent this problem we calculate its integral transform

$\begin{split}\mathrm{Im}\tilde{G}_{h}(\alpha,\beta,E)&\\\ &=\int
d\omega\mathrm{Im}G_{h}(\alpha,\beta,\omega)K(\omega,E)\\\
&=-\pi\sumint_{\Phi_{A-1}}\langle
0|a_{\beta}^{\dagger}|K\big{(}E_{\Phi},E-E_{0}\big{)}|a_{\alpha}|0\rangle\\\
&=-\pi\langle
0|a_{\beta}^{\dagger}|K\big{(}H,E-E_{0}\big{)}|a_{\alpha}|0\rangle\,,\end{split}$
(14)

with the integral kernel $K(\omega,E)$. If we are interested in
$\mathrm{Im}G_{h}(\alpha,\beta,E)$, the integral transform has to be inverted.
The inversion procedure requires in general solving an ill-posed problem. It
has been successful for the Lorentz or Laplace kernel when the responses have
a relatively simple shape composed of one or two peaks Efros  _et al._ (1994);
Raghavan  _et al._ (2021). A variety of techniques were proposed to achieve
it. The spectral function, however, has a more complicated structure, making
the inversion practically impossible. Therefore, here we follow a different
strategy outlined in Ref. Sobczyk  _et al._ (2022). We quote here only the
most important steps of the derivation, while all the details can be found in
Refs. Sobczyk  _et al._ (2022); Sobczyk and Roggero (2021).

We reconstruct $\mathrm{Im}G_{h}(\alpha,\beta,E)$ as a histogram. Our goal is
to estimate each bin (centered at $\eta$ having width $2\Delta$), as

$\mathrm{Im}G_{h}(\alpha,\beta;\eta,\Delta)\equiv\int_{\eta-\Delta}^{\eta+\Delta}dE~{}\mathrm{Im}G_{h}(\alpha,\beta,E)\,,$
(15)

using the integral transform

$\mathrm{Im}\tilde{G}_{h}(\alpha,\beta;\eta,\Delta)\equiv\int_{\eta-\Delta}^{\eta+\Delta}dE~{}\mathrm{Im}\tilde{G}_{h}(\alpha,\beta,E)\,.$
(16)

The uncertainty of this reconstruction depends on the properties of the kernel
$K$. From our previous studies, we found that the Gaussian kernel has very
convenient properties, which we characterize using parameters $\Sigma$
(accurateness) and $\Lambda$ (resolution)

$\sup_{\omega\in[-1,1]}\sumint_{\omega-\Lambda}^{\omega+\Lambda}K(\omega,E)dE\geq
1-\Sigma\,.$ (17)

Using these definitions we arrive at the histogram which is constrained from
below and above by the integrated integral transforms,

$\begin{split}\mathrm{Im}\tilde{G}_{h}(\Delta-\Lambda)-\Sigma\leq\mathrm{Im}G_{h}(\Delta)\leq\mathrm{Im}\tilde{G}_{h}(\Delta+\Lambda)+\Sigma\,,\end{split}$
(18)

where we suppressed the quantum numbers $\alpha$, $\beta$ and the bin’s center
$\eta$. The integral transform in Eq. (14) itself can be calculated in various
manners, depending on the employed kernel. For example, the Lorentz integral
transform can be conveniently obtained via Lanczos algorithm which gives
access to the set of the lowest eigenvalues. Here, however, we will use a
different strategy which can be applied not only to the Lorentzian but also to
the Gaussian kernel, namely by expanding the kernel into Chebyshev polynomials

$K(\omega,E)=\sum_{k=0}^{N}c_{k}(E)T_{k}(\omega)\,.$ (19)

We note that Chebyshev polynomials are defined on $[-1,1]$ so we have to scale
our problem in such a way that the spectrum of Hamiltonian is confined in
$[-1,1]$ range. The coefficients $c_{k}$ of this expansion have analytical
form, while the Chebyshev polynomials $T_{k}$ follow the recursive relations

$\begin{split}&T_{0}(x)=1;\,\,\,\,\,T_{-1}(x)=T_{1}(x)=x;\\\
&T_{n+1}(x)=2xT_{n}(x)-T_{n-1}(x)\,.\end{split}$ (20)

We can obtain the moments of this expansion iterating the action of the
nuclear Hamiltonian $H$ on the initial state $a_{\alpha}|0\rangle$

$\begin{split}&\langle\tilde{\Phi}_{0}|\equiv\langle
0|a_{\beta}^{\dagger}|\,,\ \ \ \ \ |\Phi_{0}\rangle\equiv
a_{\alpha}|0\rangle\,,\\\
&\langle\tilde{\Phi}_{k}|\equiv\langle\tilde{\Phi}_{k-1}|H\ \ \ \ \
|\Phi_{k}\rangle=H|\Phi_{k-1}\rangle\\\
&\mu_{0}=\langle\tilde{\Phi}_{0}|\Phi_{0}\rangle\,,\ \ \ \ \
\mu_{1}=\langle\tilde{\Phi}_{0}|\Phi_{1}\rangle\equiv\langle\tilde{\Phi}_{1}|\Phi_{0}\rangle\\\
&\mu_{k+1}=2\langle\tilde{\Phi}_{0}|\Phi_{k+1}\rangle-\mu_{k-1}\equiv
2\langle\tilde{\Phi}_{k+1}|\Phi_{0}\rangle-\mu_{k-1}\,.\end{split}$ (21)

Combining Eqs. (14) and (21) we arrive at

$\begin{split}\mathrm{Im}\tilde{G}_{h}(\alpha,\beta,E)=&-\pi\sum_{k=0}^{N}c_{k}(E)\langle
0|a_{\beta}^{\dagger}T_{k}\big{(}H\big{)}a_{\alpha}|0\rangle\\\
\equiv&-\pi\sum_{k=0}^{N}c_{k}(E)\mu_{k}\,.\end{split}$ (22)

We truncate the expansion at the level on $N$ moments, introducing a
controllable error $\gamma$,

$\gamma=\sum_{k=N+1}^{\infty}c_{k}(E)T_{k}(\omega)\,,$ (23)

It has to be included into an overall uncertainty budget, leading to the final
prediction

$\begin{split}\mathrm{Im}\tilde{G}_{h}(\Delta-&\Lambda)-\Sigma-2\gamma(\Delta-\Lambda)\\\
&\leq\mathrm{Im}G_{h}(\Delta)\\\
\leq\mathrm{Im}\tilde{G}_{h}&(\Delta+\Lambda)+\Sigma+2\gamma(\Delta+\Lambda)\,.\end{split}$
(24)

By setting $\Delta$ (the histogram’s width), $\Lambda$ (width of the kernel)
and $N$ (number of Chebyshev moments) we can estimate the lower and upper
bound on $\mathrm{Im}G_{h}(\Delta)$ according to Eq. (24). We note that
$\Sigma$ and $\gamma$ have a known analytical form, and the uncertainty is
mainly driven by
$|\mathrm{Im}\tilde{G}_{h}(\Delta+\Lambda)-\mathrm{Im}\tilde{G}_{h}(\Delta-\Lambda)|$.

### Coupled-cluster theory

We calculate the Chebyshev moments $\mu_{k}$ within the spherical coupled-
cluster framework Hagen  _et al._ (2014). This formalism starts from a
reference state $|\Phi\rangle$, in our case a Hartree-Fock solution, on top of
which we include the nuclear correlations using an exponential ansatz

$\begin{split}&|0\rangle=e^{T}|\Phi\rangle\,,\\\
&T=\sum_{i,a}t^{a}_{i}a^{\dagger}_{a}a_{i}+\frac{1}{4}\sum_{ijab}t^{ab}_{ij}a_{a}^{\dagger}a_{b}^{\dagger}a_{i}a_{j}+...\,.\end{split}$
(25)

In our present calculation we retain the first two terms of this expansion,
i.e. we work in the singles and doubles (CCSD) approximation. The $t$
amplitudes appearing in the correlation operator $T$ can be determined solving
a set of coupled nonlinear equations. For the calculation of the Green’s
function we need to construct a set of initial $A-1$ states acting with the
similarity transformed operators $\overline{a}_{\alpha}$ and
$\overline{a}^{\dagger}_{\alpha}$ on the left and right ground-state

$\begin{split}&|\Phi_{0}\rangle=\overline{a}_{\alpha}|\Phi\rangle\equiv
e^{-T}a_{\alpha}e^{T}|\Phi\rangle\,,\\\
&\langle\tilde{\Phi}_{0}|=\langle\Phi|\overline{a}^{\dagger}_{\alpha}\equiv\langle\Phi|(1+\Lambda)e^{-T}a^{\dagger}_{\alpha}e^{T}\,,\end{split}$
(26)

where $\Lambda$ is the de-excitation operator which has to be included since
the coupled-cluster is a non-hermitian theory yielding different left and
right eigenstates. The calculation of Chebyshev moments requires an iterative
action of the similarity transformed Hamiltonian, following the recursive
relation from Eq. (21).

##  IV Results

In all the results presented in this paper, we employ the NNLOsat nuclear
Hamiltonian Ekström  _et al._ (2015) containing both nucleon-nucleon (NN) and
three-nucleon (3N) forces derived in chiral effective field theory at next-to-
next-to leading order Epelbaum  _et al._ (2009). The low-energy constants in
this Hamiltonian are fitted both to NN scattering data, as well as to
properties of light nuclei and selected medium-mass nuclei. We recall that the
current operators implemented in this work are not derived in chiral effective
field theory, but we rather use the relativistic forms described in Section
II. In the coupled-cluster calculations, 3N interactions are approximated at
the normal-ordered two-body level Hagen  _et al._ (2007); Roth  _et al._
(2012), and an additional cut on three-nucleon configurations $E_{3max}\leq
16$ is imposed. We performed calculations for the model space of 15 oscillator
shells and values of underlying harmonic oscillator frequencies
$\hbar\Omega=12-20$ MeV.

To benchmark our calculation we first look at the charge distribution and
compare it with previous results from the SCGF Rocco and Barbieri (2018)
method. In the upper panel of Fig. 2 we present the direct result of the
computation which includes spurious center of mass (CoM) contaminations,
denoted with CCSD. We also show the intrinsic charge distribution, for which
the CoM contributions were subtracted (see Ref. Sobczyk  _et al._ (2022) for
details), denoted with “CCSD intr”. They are both in a very good agreement
with the SCGF predictions, for which a different numerical procedure is used
to remove the CoM contributions. We also note that the “CCSD” and “CCSD intr”
distributions are similar, which confirms the well known fact that spurious
CoM effects decrease with the nuclear mass (for a comparison see Fig. 2 in
Ref. Sobczyk  _et al._ (2022) for 4He, where the effect was larger).

Next, we look at the momentum distribution $n(\mathbf{p})=\int
dE\,S(\mathbf{p},E)$ to further assess the role of the CoM contamination We
follow the same procedure as explained in Ref. Sobczyk  _et al._ (2022) to
calculate the intrinsic $n(\mathbf{p})$. In Fig. 2 we show both
$n(\mathbf{p})$ and $\mathbf{p}^{2}n(\mathbf{p})$ in the inset. The CoM
removal affects mostly low momenta (shown as a blue band in Fig. 2). The
uncertainty comes from varying the width of the CoM Gaussian used as an
ansatz, $\hbar\tilde{\Omega}=16-24$ MeV. The CoM effect will be negligible
when we consider the momentum-weighted $\mathbf{p}^{2}n(\mathbf{p})$, as in
the case of cross-section calculation. We will therefore safely neglect this
effect from now on.

Figure 2: Charge distribution (upper panel) comparison between CCSD and SCGF
calculations using the same interaction NNLOsat. Momentum distribution (lower
panel) using CCSD. See text for details.

In order to investigate the dependence of the spectral function on the basis,
we looked separately at integrated distributions: $n(\mathbf{p})=\int dE\
S(\mathbf{p},E)$ and $S(E)=\int d^{3}\mathbf{p}\ S(\mathbf{p},E)$. Momentum
distribution $n(\mathbf{p})$ is practically independent on the choice of
$\hbar\Omega$. The differences between various $S(E)$ can be appreciated in
Fig. 3. The dominating peaks (below $30$ MeV) have almost the same strength
and are shifted by less than $1$ MeV, while the spectrum above 30 MeV is quite
different, as can be seen in the inset of Fig. 3. A direct comparison of the
full 2D distribution of the spectral function for various $\hbar\Omega$
reveals some strength redistribution. There are also some small regions which
give a negative contribution. This non-physical behaviour is most likely due
to the fact that the coupled-cluster theory is non-hermitian. The appearance
of small admixtures of non-physical states has been already observed Gu  _et
al._ (2023); Hagen  _et al._ (2010). We have numerically checked that the
negative contribution is smallest for $\hbar\Omega=14$ MeV and in this case it
stays at the per-mil level. For other values of oscillator frequencies it
reaches at most $4\%$ for $\hbar\Omega=20$ MeV. Therefore, we have decided to
perform all the further calculations with this optimal value of
$\hbar\Omega=14$ MeV which alleviates the non-physical behaviour.

Figure 3: The integrated energy distribution $S(E)$ of the proton spectral
function for various values of $\hbar\Omega=12-20$ MeV. The chosen binning is
also shown (arbitrary normalization).

In Fig. 4, we show the final SF histograms separately for protons and neutrons
using $2\Delta=3.5$ MeV binning. In this case the Hamiltonian spectrum is
limited by energy $E=250$ MeV.222This value is needed to scale the spectrum to
$[-1,1]$ range where the Chebyshev polynomials are defined. We set
$\Lambda=0.66$ MeV, and the number of Chebyshev moments $N=4000$ to keep the
truncation error $\gamma$ negligibly small. One can observe two clearly
dominating peaks at $E\approx 13(10)$ and $21(18)$ MeV for neutrons (protons),
which correspond to $1p_{1/2}$ and $1p_{3/2}$ states, and some strength
distributed at higher energies. We note that for this estimation we use
$\mathrm{Im}G_{h}(\Delta)\approx\mathrm{Im}\tilde{G}_{h}(\Delta)$ which lies
between the lower and upper bounds, as shown in Eq. (24).

Figure 4: 16O spectral functions for protons (upper panel) and neutrons (lower
panel). See text for details of spectral reconstruction. The theoretical
uncertainty, not shown in this figure, is discussed in Subsec. IV.1.

###  IV.1 Uncertainty estimation

The estimated uncertainties coming from the ChEK procedure (see Eq. (24))
affect mostly the energy range $E>40$ MeV. This can be understood when various
sources of uncertainty are analyzed in Eq. (24). With our choice of parameters
we keep $\Sigma$ and $\gamma$ small, and the uncertainty is driven by
$|\mathrm{Im}\tilde{G}_{h}(\Delta+\Lambda)-\mathrm{Im}\tilde{G}_{h}(\Delta-\Lambda)|$.
The lower part of the spectrum (below 30 MeV) is composed of well separated
peaks, and therefore with our choice of the histogram binning the
uncertainties are negligible. The total strength of the SF is dominated by
this region. Therefore, the uncertainties have an overall small impact on the
cross section. The uncertainties estimated for
$\mathrm{Im}\tilde{G}(\alpha,\beta,E)$ lead to $|\mathbf{p}|$-dependent errors
in the final SF, according to Eq. (11). In Fig. 5 we show the proton spectral
function at three values of momenta $|\mathbf{p}|$. The spectrum below $E=40$
MeV is not affected by the uncertainties while for higher energies the errors
are larger. They will not, however, influence much the cross-section results,
since the hadron tensor is weighted by $\mathbf{p}^{2}d|\mathbf{p}|$. In fact,
the weighted momentum distribution $\mathbf{p}^{2}n(\mathbf{p})$ peaks at
$|\mathbf{p}|\approx 150$ MeV (see Fig. 2).

Figure 5: Proton spectral function of 16O for three values of momentum
$|\mathbf{p}|$. The uncertainty bars come from Eq. (24).

###  IV.2 Applications to electron-nucleus scattering

Electron scattering experiments serve as an excellent test to check the
reliability of the nuclear models and of the assumed approximations.
Unfortunately, there are only scarce data available for 16O, corresponding to
momentum transfers in the range of $320-650$ MeV.

Figure 6: Electron scattering differential cross-section on 16O for different
kinematics which correspond to the momentum transfer $|\mathbf{q}|\approx
320-650$ MeV. We show results obtained with CCSD without including the optical
potential (dashed line) and after its inclusion (solid line). The uncertainty
bands come from the uncertainty of the ChEK method. Experimental data was
taken from Ref. O’Connell  _et al._ (1987); Anghinolfi  _et al._ (1996).

In Fig. 6, we compare our results with the available experimental data.
Clearly, our calculations compare better in the kinematical regime of higher
momentum transfer, as expected from an IA assumption. When accounting for the
FSI via the optical potential as described in Sec. II, the position of the QE
peak is shifted to lower energy transfers so that the agreement with the
experimental data improves. As anticipated, the effect of FSI is stronger in
the lower momentum-energy regimes presented in the first row of Fig. 6.

Within our approach, we propagate the theoretical uncertainty from the SF to
the cross-section results. In practice, we construct two SFs taking both the
lowest and the highest values for each histogram bin as presented in Fig. 5,
and with those two SFs we construct a lower and an upper cross section,
respectively, which lead to the bands in Fig. 6. The obtained uncertainty is
of the order of a few percents and reaches about the $10\%$ mark at the QE
peak. Although the response in the QE peak for $|\mathbf{q}|>500$ MeV seems to
well describe the data, our predictions might actually be too high. In fact,
other mechanisms not included in the calculation, such as meson exchange
currents (MEC) and pion production, typically contribute by mostly enhancing
the high-energy slope of the QE peak. However, we note that the imaginary part
of the optical potential, not included in our current calculations, could
quench the response, giving therefore some room for the above mentioned
enhancing contributions. Therefore, we expect here some cancellations of the
omitted contributions, whose investigation is left to future work.

Finally, in Fig. 7 we present a comparison of our calculations with the
previous results of Ref. Rocco and Barbieri (2018) obtained from the SCGF
method. We show only one kinematics, since the comparison is similar for other
setups. We observe that the two curves are very similar, indicating a nice
agreement. Looking at the details, the QE peak in CCSD is slightly shifted
towards the left with respect to the SCGF calculation. We expect that the
source of this tiny deviation lies in the differences between the two many-
body methods. As pointed out before, we obtain also very similar charge
distributions. We have checked that the momentum distributions are also in
very good agreement between the two methods, indicating that the many-body
description of the ground states is very similar. However, the slight
differences observed in the cross section may indicate a stronger sensitivity
to the details of the many-body method for dynamical observables. We
nevertheless consider this benchmark very successful.

Figure 7: Electron scattering cross section for 16O calculated within CCSD
(continuous blue line) and within SCGF Rocco and Barbieri (2018) (black
dashed-dotted line) in comparison to experimental data from Ref. Anghinolfi
_et al._ (1996).

###  IV.3 Applications to neutrino-nucleus scattering

We now redirect our attention to the application of our calculations to
neutrino-nucleus scattering. The 16O spectral function is of particular
interest for T2K and future T2HK experiments. A direct comparison of the QE
peak with the data for the neutrino-nucleus scattering is currently not
possible due to the experimental constraints. The neutrino flux has a broad
energy spectrum so that many mechanisms contribute and cannot be well
separated. Moreover, the current experimental uncertainties are large and
dominated by the statistics. Nevertheless, recently the T2K collaboration
published inclusive cross-section on 16O for CC$0\pi$ events (no pions
detected in the final state) for various angles of the outgoing muon Abe  _et
al._ (2020). This observable should have a large contribution coming from the
single-nucleon knockout mechanism, which is well described by our SF. To make
a full comparison with the data, we need to $(i)$ account for all other
possible mechanisms, beyond the one-nucleon knockout and $(ii)$ get the
distribution of the produced hadrons going beyond the inclusive cross-section.
We achieve this by implementing our SF in the NuWro MC event generator
Juszczak  _et al._ (2006); Golan  _et al._ (2012). Typically, the MC
generators describe the neutrino-nucleus scattering in a two-step process. In
the first step, the scattering takes place on a single nucleon (or a pair of
nucleons in case of meson-exchange currents) in the primary vertex. This is
where we include our SF model. In the next step, the produced particles
(predominantly nucleons and pions) are cascaded through the nucleus, where
they can re-scatter, be absorbed, or produce other particles. Various
approaches to model the inter-nuclear cascade were recently compared in Ref.
Dytman  _et al._ (2021). From this perspective including the imaginary part of
the optical potential – omitted in our calculation – might lead to the double
counting of some effects already accounted for in the cascade.

The results of the simulation for the double differential cross section
$\nu_{\mu}+^{16}\mathrm{O}\to\mu^{-}+X$ of CC0$\pi$ events done with NuWro and
our SF are presented in Fig. 8. In the same plot we also show the predictions
for the SF with the inclusion of optical potential (hatched pattern, denoted
with SF+FSI). Here, we do not account for the theoretical uncertainty of
spectral function, focusing only on the role played by optical potential.
Other dynamical channels, MEC and resonance contributions (RES), were chosen
to be the same as explained in Ref. Abe  _et al._ (2020). We are aware that
these predictions are not fully consistent, since the theoretical description
of each mechanism is based on a different model. To make our comparison with
the experimental data more meaningful, we would need to address all the
contributions within the same SF method. This is certainly an important
direction of future investigations. However, at this point we focus only on
the IA mechanism. In Fig. 9 we show our final prediction, including the
uncertainty of SF. It leads up to $\sim 10\%$ effect, depending on the
considered kinematics. We find a reasonable agreement with the data, similar
to the results of Ref. Abe  _et al._ (2020) where several MC event generators
were employed (a variety of models were used for the QE mechanism, including
random phase approximation corrections, a phenomenological spectral function
or relativistic mean field calculations). Our SF quenches the response when
compared to the local Fermi Gas, even by $25\%$ for forward scattering angles.
In fact, the kinematics of the most forward angles (upper left panel in Fig.
8) depend mostly on the details of the nuclear model, since the momentum
transfer is the smallest (covers mostly the range of $|\mathbf{q}|\approx
100-300$ MeV). For the angles $0.96<\cos\theta<1$ the optical potential causes
a substantial depletion of the bins corresponding to the values of
$|\mathbf{k}^{\prime}|<800$ MeV. Our simple model of FSI gives reasonable
results for the electron scattering at the intermediate momentum transfer
$|\mathbf{q}|\gtrapprox 450$ MeV shifting the QE peak to lower energies (see
Fig. 6). Here, on the contrary, we observe a significant effect. At this
kinematics (low $|\mathbf{q}|$ and forward scattering angles) the IA is much
less reliable. This region of phase-space escapes the capability of our method
and should be rather described by a consistent calculation of FSI, available
with the LIT-CC approach. We also observe that the contribution coming from
the phenomenological MEC is substantial in this range. This prediction can
also be verified using an ab-initio approach including one- and two-body
currents Lovato  _et al._ (2020); Pastore  _et al._ (2020).

Figure 8: The double differential cross section
$\nu_{\mu}+^{16}\mathrm{O}\to\mu^{-}+X$ of CC0$\pi$ events measured by the T2K
experiment Abe  _et al._ (2020). The momentum distribution
$|\mathbf{k}^{\prime}|$ of outgoing $\mu^{-}$ was measured for five ranges of
the scattering angle $\cos\theta$. The results “SF” were obtained using our
spectral function, while “SF+FSI” include also the real part of the optical
potential. In both cases we do not include theoretical uncertainties. Other
mechanisms, predominantly MEC and resonance production (RES) give smaller
contribution according to the model used in the simulation Abe  _et al._
(2020).

Figure 9: The same as Fig. 8, showing theoretical uncertainty coming from the
reconstruction of spectral function. “SF+FSI (low)” was obtained using the
lower bound of the spectral reconstruction, while “SF+FSI (high)” corresponds
to the upper bound.

##  V Conclusion and outlook

We calculated spectral functions of 16O within the many-body coupled-cluster
framework and employing a chiral nuclear Hamiltonian including 3N forces at
next-to-next-to leading order. The procedure required a reconstruction of the
spectral properties (i.e. the energy-dependant part of the SF), which we
performed within the ChEK method. This approach, which was benchmarked on the
4He in an earlier publication Sobczyk  _et al._ (2022), allows to assess the
uncertainty of our calculation and to propagate it to the cross-section
results.

Within the impulse approximation, the SF can be directly related to the
scattering cross-section. Using this assumption, we give predictions for the
lepton-nucleus scattering in the QE regime both for electron and neutrino
scattering. The electron scattering data for 16O are scarce and cover only the
medium and high momentum transfer regions.

Within this energy range the IA works well and we get a good agreement with
the data, although for $|\mathbf{q}|<500$ MeV the FSI play an important role
and the inclusion of optical potential visibly improves the agreement. We
still do not account for the absorption of the outgoing nucleon, which is
certainly an important topic to be addressed in the future when aiming at the
comparison with inclusive data. We would like to point out that further
investigations of the QE region are currently restricted due to the lack of
low-energy electron-scattering data on 16O. More data would be of great value
to guide theoretical models used in the future T2HK experiment. There are
plans to take new data on 16O in the future at MAMI in Germany Ankowski  _et
al._ (2022).

We presented a comparison with neutrino T2K data for CC$0\pi$ events which are
sensitive to the QE mechanism. To this end, we implemented our SF in the NuWro
MC generator. In our analysis we observed that for the forward angles the
optical potential plays an important role. In fact, the IA picture becomes
much less reliable in this regime and a consistent calculation which accounts
for the final state interactions, as the LIT-CC, would be more appropriate.
Also the role played by the two-body currents should be examined. The work in
this direction is already on-going. Since our studies are mainly motivated by
the neutrino oscillation experiments, we find it important to make our
spectral functions available for further exploration rep .

###### Acknowledgements.

We acknowledge useful discussions with G. Hagen and T. Papenbrock and we thank
them for letting us use the NuCCore coupled-cluster code. J.E.S. acknowledges
the support of the Humboldt Foundation through a Humboldt Research Fellowship
for Postdoctoral Researchers. This project has received funding from the
European Union’s Horizon 2020 research and innovation programme under the
Marie Skłodowska-Curie grant agreement No. 101026014. This work was supported
in part by the Deutsche Forschungsgemeinschaft (DFG) through the Cluster of
Excellence “Precision Physics, Fundamental Interactions, and Structure of
Matter” (PRISMA+ EXC 2118/1) funded by the DFG within the German Excellence
Strategy (Project ID 39083149)

## References

  * Acciarri  _et al._ (2015) R. Acciarri _et al._ (DUNE), “Long-Baseline Neutrino Facility (LBNF) and Deep Underground Neutrino Experiment (DUNE),” (2015), arXiv:1512.06148 [physics.ins-det] .
  * Abe  _et al._ (2015) K. Abe _et al._ (Hyper-Kamiokande Proto-Collaboration), “Physics potential of a long-baseline neutrino oscillation experiment using a J-PARC neutrino beam and Hyper-Kamiokande,” PTEP 2015, 053C02 (2015).
  * Alvarez-Ruso  _et al._ (2018) L. Alvarez-Ruso _et al._ (NuSTEC), “NuSTEC White Paper: Status and challenges of neutrino–nucleus scattering,” Prog. Part. Nucl. Phys. 100, 1–68 (2018), arXiv:1706.03621 [hep-ph] .
  * Ankowski  _et al._ (2022) A. M. Ankowski _et al._ , “Electron Scattering and Neutrino Physics,” (2022), arXiv:2203.06853 [hep-ex] .
  * Ruso  _et al._ (2022) L. Alvarez Ruso _et al._ , “Theoretical tools for neutrino scattering: interplay between lattice QCD, EFTs, nuclear physics, phenomenology, and neutrino event generators,” (2022), arXiv:2203.09030 [hep-ph] .
  * Lovato  _et al._ (2018) A. Lovato, S. Gandolfi, J. Carlson, Ewing Lusk, Steven C. Pieper, and R. Schiavilla, “Quantum Monte Carlo calculation of neutral-current $\nu-^{12}C$ inclusive quasielastic scattering,” Phys. Rev. C 97, 022502 (2018), arXiv:1711.02047 [nucl-th] .
  * Lovato  _et al._ (2020) A. Lovato, J. Carlson, S. Gandolfi, N. Rocco, and R. Schiavilla, “Ab initio study of $({\nu}_{\ell},{\ell}^{-})$ and $({\overline{\nu}}_{\ell},{\ell}^{+})$ inclusive scattering in ${}^{12}\mathrm{C}$: Confronting the miniboone and t2k ccqe data,” Phys. Rev. X 10, 031068 (2020).
  * Sobczyk  _et al._ (2020) J. E. Sobczyk, B. Acharya, S. Bacca, and G. Hagen, “Coulomb sum rule for 4He and 16O from coupled-cluster theory,” Phys. Rev. C 102, 064312 (2020), arXiv:2009.01761 [nucl-th] .
  * Sobczyk  _et al._ (2021) J. E. Sobczyk, B. Acharya, S. Bacca, and G. Hagen, “Ab initio computation of the longitudinal response function in 40Ca,” Phys. Rev. Lett. 127, 072501 (2021), arXiv:2103.06786 [nucl-th] .
  * Rocco  _et al._ (2019) Noemi Rocco, Satoshi X. Nakamura, T. S. H. Lee, and Alessandro Lovato, “Electroweak Pion-Production on Nuclei within the Extended Factorization Scheme,” Phys. Rev. C 100, 045503 (2019), arXiv:1907.01093 [nucl-th] .
  * Jiang  _et al._ (2022) L. Jiang _et al._ (Jefferson Lab Hall A), “Determination of the argon spectral function from (e,e’p) data,” Phys. Rev. D 105, 112002 (2022), arXiv:2203.01748 [nucl-ex] .
  * Benhar  _et al._ (1994) O. Benhar, A. Fabrocini, S. Fantoni, and I. Sick, “Spectral function of finite nuclei and scattering of GeV electrons,” Nucl. Phys. A 579, 493–517 (1994).
  * Nieves and Sobczyk (2017) Juan Nieves and Joanna Ewa Sobczyk, “In medium dispersion relation effects in nuclear inclusive reactions at intermediate and low energies,” Annals Phys. 383, 455–496 (2017), arXiv:1701.03628 [nucl-th] .
  * Buss  _et al._ (2007) O. Buss, T. Leitner, U. Mosel, and L. Alvarez-Ruso, “The Influence of the nuclear medium on inclusive electron and neutrino scattering off nuclei,” Phys. Rev. C 76, 035502 (2007), arXiv:0707.0232 [nucl-th] .
  * Rocco and Barbieri (2018) N. Rocco and C. Barbieri, “Inclusive electron-nucleus cross section within the Self Consistent Green’s Function approach,” Phys. Rev. C 98, 025501 (2018), arXiv:1803.00825 [nucl-th] .
  * Barbieri  _et al._ (2019) C. Barbieri, N. Rocco, and V. Somà, “Lepton Scattering from 40Ar and Ti in the Quasielastic Peak Region,” Phys. Rev. C 100, 062501 (2019), arXiv:1907.01122 [nucl-th] .
  * Roggero (2020) Alessandro Roggero, “Spectral density estimation with the Gaussian Integral Transform,” Phys. Rev. A 102, 022409 (2020), arXiv:2004.04889 [quant-ph] .
  * Sobczyk and Roggero (2021) Joanna E. Sobczyk and Alessandro Roggero, “Spectral density reconstruction with Chebyshev polynomials,” (2021), arXiv:2110.02108 [nucl-th] .
  * Hagen  _et al._ (2014) G. Hagen, T. Papenbrock, M. Hjorth-Jensen, and D. J. Dean, “Coupled-cluster computations of atomic nuclei,” Rep. Prog. Phys. 77, 096302 (2014).
  * Sobczyk  _et al._ (2022) J. E. Sobczyk, S. Bacca, G. Hagen, and T. Papenbrock, “Spectral function for 4He using the Chebyshev expansion in coupled-cluster theory,” (2022), arXiv:2205.03592 [nucl-th] .
  * Juszczak  _et al._ (2006) Cezary Juszczak, Jaroslaw A. Nowak, and Jan T. Sobczyk, “Simulations from a new neutrino event generator,” Nucl. Phys. B Proc. Suppl. 159, 211–216 (2006), arXiv:hep-ph/0512365 .
  * Golan  _et al._ (2012) T. Golan, J. T. Sobczyk, and J. Zmuda, “NuWro: the Wroclaw Monte Carlo Generator of Neutrino Interactions,” Nucl. Phys. B Proc. Suppl. 229-232, 499–499 (2012).
  * Bradford  _et al._ (2006) R. Bradford, A. Bodek, Howard Scott Budd, and J. Arrington, “A New parameterization of the nucleon elastic form-factors,” Nucl. Phys. B Proc. Suppl. 159, 127–132 (2006), arXiv:hep-ex/0602017 .
  * Sobczyk (2017) Joanna Ewa Sobczyk, “Intercomparison of lepton-nucleus scattering models in the quasielastic region,” Phys. Rev. C 96, 045501 (2017), arXiv:1706.06739 [nucl-th] .
  * Cooper  _et al._ (1993) E. D. Cooper, S. Hama, B. C. Clark, and R. L. Mercer, “Global Dirac phenomenology for proton nucleus elastic scattering,” Phys. Rev. C 47, 297–311 (1993).
  * Ankowski  _et al._ (2015) Artur M. Ankowski, Omar Benhar, and Makoto Sakuda, “Improving the accuracy of neutrino energy reconstruction in charged-current quasielastic scattering off nuclear targets,” Phys. Rev. D 91, 033005 (2015), arXiv:1404.5687 [nucl-th] .
  * Benhar  _et al._ (1991) O. Benhar, A. Fabrocini, S. Fantoni, G. A. Miller, V. R. Pandharipande, and I. Sick, “Scattering of GeV electrons by nuclear matter,” Phys. Rev. C 44, 2328–2342 (1991).
  * Efros  _et al._ (1994) Victor D. Efros, Winfried Leidemann, and Giuseppina Orlandini, “Response functions from integral transforms with a Lorentz kernel,” Phys. Lett. B 338, 130–133 (1994), arXiv:nucl-th/9409004 .
  * Raghavan  _et al._ (2021) Krishnan Raghavan, Prasanna Balaprakash, Alessandro Lovato, Noemi Rocco, and Stefan M. Wild, “Machine learning-based inversion of nuclear responses,” Phys. Rev. C 103, 035502 (2021), arXiv:2010.12703 [nucl-th] .
  * Ekström  _et al._ (2015) A. Ekström, G. R. Jansen, K. A. Wendt, G. Hagen, T. Papenbrock, B. D. Carlsson, C. Forssén, M. Hjorth-Jensen, P. Navrátil, and W. Nazarewicz, “Accurate nuclear radii and binding energies from a chiral interaction,” Phys. Rev. C 91, 051301 (2015), arXiv:1502.04682 [nucl-th] .
  * Epelbaum  _et al._ (2009) Evgeny Epelbaum, Hans-Werner Hammer, and Ulf-G. Meissner, “Modern Theory of Nuclear Forces,” Rev. Mod. Phys. 81, 1773–1825 (2009), arXiv:0811.1338 [nucl-th] .
  * Hagen  _et al._ (2007) G. Hagen, T. Papenbrock, D. J. Dean, A. Schwenk, A. Nogga, M. Wloch, and P. Piecuch, “Coupled-cluster theory for three-body Hamiltonians,” Phys. Rev. C 76, 034302 (2007), arXiv:0704.2854 [nucl-th] .
  * Roth  _et al._ (2012) Robert Roth, Sven Binder, Klaus Vobig, Angelo Calci, Joachim Langhammer, and Petr Navratil, “Ab Initio Calculations of Medium-Mass Nuclei with Normal-Ordered Chiral NN+3N Interactions,” Phys. Rev. Lett. 109, 052501 (2012), arXiv:1112.0287 [nucl-th] .
  * Gu  _et al._ (2023) Chenyi Gu, Z. H. Sun, G. Hagen, and T. Papenbrock, “Entanglement entropy of nuclear systems,” (2023), arXiv:2303.04799 [nucl-th] .
  * Hagen  _et al._ (2010) G. Hagen, T. Papenbrock, D. J. Dean, and M. Hjorth-Jensen, “Ab initio coupled-cluster approach to nuclear structure with modern nucleon-nucleon interactions,” Phys. Rev. C 82, 034330 (2010), arXiv:1005.2627 [nucl-th] .
  * O’Connell  _et al._ (1987) J. S. O’Connell _et al._ , “Electromagnetic excitation of the delta resonance in nuclei,” Phys. Rev. C35, 1063 (1987).
  * Anghinolfi  _et al._ (1996) M. Anghinolfi _et al._ , “Quasi-elastic and inelastic inclusive electron scattering from an oxygen jet target,” Nucl. Phys. A602, 405–422 (1996), nucl-th/9603001 .
  * Abe  _et al._ (2020) K. Abe _et al._ (T2K), “Simultaneous measurement of the muon neutrino charged-current cross section on oxygen and carbon without pions in the final state at T2K,” Phys. Rev. D 101, 112004 (2020), arXiv:2004.05434 [hep-ex] .
  * Dytman  _et al._ (2021) Steven Dytman, Yoshinari Hayato, Roland Raboanary, J. T. Sobczyk, Julia Tena Vidal, and Narisoa Vololoniaina, “Comparison of validation methods of simulations for final state interactions in hadron production experiments,” Phys. Rev. D 104, 053006 (2021), arXiv:2103.07535 [hep-ph] .
  * Pastore  _et al._ (2020) Saori Pastore, Joseph Carlson, Stefano Gandolfi, Rocco Schiavilla, and Robert B. Wiringa, “Quasielastic lepton scattering and back-to-back nucleons in the short-time approximation,” Phys. Rev. C 101, 044612 (2020), arXiv:1909.06400 [nucl-th] .
  * (41) https://github.com/jesobczyk/spectral-function.