# Point Spread Function Deconvolution Using a Convolutional Autoencoder for
Astronomical Applications

Sreevarsha Sreejith1,3 Anže Slosar1 Hong Wang2 1 Physics Department,
Brookhaven National Laboratory, Upton, NY 11973  2 Computational Science
Initiative, Brookhaven National Laboratory, Upton, NY 11973  3 Physics
Department, University of Surrey, Guildford GU2 7XH, UK

###### Abstract

A major issue in optical astronomical image analysis is the combined effect of
the instrument’s point spread function (PSF) and the atmospheric seeing that
blurs images and changes their shape in a way that is band and time-of-
observation dependent. In this work we present a very simple neural network
based approach to non-blind image deconvolution that relies on feeding a
Convolutional Autoencoder (CAE) input images that have been preprocessed by
convolution with the corresponding PSF and its regularized inverse. Compared
to our previous work based on Deep Wiener Deconvolution, the new approach is
conceptually simpler and computationally much less intensive while achieving
only marginally worse results. In this work we also present a new approach for
dealing with limited input dynamic range of neural networks compared to the
dynamic range present in astronomical images.

##  1 Introduction

Astronomy is entering into an era of big data with existing and upcoming
ground-based and space-borne observatories probing remarkably large volumes of
the observable sky with depths and cadences that are hitherto unseen. With
upcoming surveys such as the Legacy Survey of Space and Time (LSST, LSST
Science Collaboration et al. 2009) to be conducted at the Vera C. Rubin
observatory, combined with existing data sets from the Dark Energy Survey
(DES, Abbott et al. 2018) & Hyper-Suprime Cam (HSC, Hikage et al. 2019), and
space-based surveys to be conducted by James Webb (JWST, Gardner et al. 2006),
Euclid (Refregier et al. 2010) & Nancy Grace Roman (Green et al. 2012)
telescopes, this ‘big data overload’ is only expected to become more
overwhelming. Both space-borne and ground-based telescopes have certain
hurdles to overcome in terms of resolution. While space-based observatories
are diffraction limited, ground-based telescopes have to cope with atmospheric
seeing. Atmospheric seeing occurs when an optical wavefront passes through
atmospheric turbulence and the perturbations in the wavefront causes the image
to get distorted, resulting in a finite point spread function (PSF). Since the
field of view is large, it is not yet feasible to implement adaptive optics
across the focal plane, leading to a PSF that is arcsecond-sized.

Since the atmosphere surrounding the telescope is variable, the PSF is
expected to vary in size, shape and orientation from one observation to
another. It also varies according to the passbands used for imaging and is
generally asymmetric. Fortunately, PSF is precisely known for each exposure
because stars contained in each exposure are excellent point sources and can
be used to directly measure the PSF.

In the face of varying PSF and noise, we want to recover unbiased shapes,
positions and fluxes of astronomical objects and machine learning techniques
have been used in recent years to great effect for this. Most of these methods
have been focused on galaxy deblending. Reiman & Göhre (2019) proposed a
branched deblender that makes use of generative adversarial networks (GANs) to
deblend overlapping field galaxies. Boucaud et al. (2020) explored the usage
of a simple convolutional neural network (CNN) in conjunction with a U-Net
(Ronneberger et al. 2015) in order to do image segmentation/galaxy deblending
and measure photometry. Arcelin et al. (2021) developed a variational
autoencoder (VAE) - like network and we ourselves explored using residual
dense neural networks (RDN) for the same purpose (Wang et al. 2022).

Most of the existing approaches are formulated with the assumption of a
constant PSF, with the caveat that the model could be retrained/modified with
techniques like transfer learning with different data sets with varying PSFs.
This could be impractical for the reasons that we have laid out in our earlier
paper, Wang et al. 2023 (HW23 henceforth). Moreover, the assumption of a
constant PSF could lead to biased photometry for individual objects. Therefore
in order to apply neural network approaches to real world scenarios, it is
necessary to address the issue of how to treat image PSFs (in addition to
other problems such as artefacts, blending, masking etc.). In this work, as in
HW23, we focus specifically on tackling the PSF issue and propose a simple
neural network architecture that, with minimal image preprocessing, removes
the residual PSF dependence and recovers object positions and shapes
effectively.

##  2 Synthetic Data

The base data set that is used in this work was simulated in the same manner
as in HW23, using the deep generative models proposed in Lanusse et al.
(2021). These models are generated using a combination of a hybrid variational
autoencoder (Kingma & Welling 2013) with the aggregate posterior distribution
modelled by a latent-space normalizing flow, termed as Flow-VAE, and was
trained on a data set based on the HST/ACS COSMOS survey (Koekemoer et al.
2007; Scoville et al. 2007) and rendered using GalSim (Rowe et al. 2015). The
physical parameters that are to be specified for image generation are half-
light radius (flux_radius, as a proxy for size), apparent magnitude in the
i-band (mag_auto, derived from Source Extractor, Bertin & Arnouts 1996) and
photometric redshift (zphot). The values for these parameters were drawn from
the respective uniform distributions as, $5\leq\texttt{flux\\_radius}\leq 15$,
$5\leq\texttt{mag\\_auto}\leq 25$ and $0\leq\texttt{zphot}\leq 2$. In order to
make the galaxy placement in the images more realistic, we slightly offset the
objects from the centre in both x and y directions with a value randomly drawn
from a uniform distribution between (-5,5).

The first major difference from the data set used in HW23 is that we do not
convolve the generated galaxy image with a small non-functional PSF that
aliases and removes the modes that might lead to over-deconvolution. Since the
deep generative model only provides a single band (band 1) and we are
interested in multi band deconvolutions, we generate two more bands with
simple non-linear transformations of band 1 (see Equation 24 and the $\alpha$
and $\beta$ values in HW23) and pile them along the third axis to form a
3-band image, which makes up the ‘truth’ data set for this work. Then each
individual band of the truth data set is convolved with a random Moffat PSF
($0.6\leq\texttt{FWHM}\leq 1$ arcsecond, $2\leq\beta\leq 5$, $-0.9\leq
g1,g2\leq 0.9$, $\sigma_{g1,g2}=0.1$, $\mu_{g1,g2}=0$) and then random
Gaussian noise ($\texttt{mean}=0$, $0\leq\texttt{variance}\leq 0.1$,
$\sigma=\texttt{variance}^{0.25}$ ) was added to create the noisy images. This
is the second major difference from HW23 –instead of using the addNoiseSNR
function from Galsim, we manually add noise to the images. The pixel scale of
the PSF of band 1 is set as the pixel scale of the 3-band image which varies
for each object in the data set, with a range between 0.088 and 0.255
pixels/arcseconds. The signal-to-noise ratio (SNR) of the data set is
calculated as Equation 25 in HW23 and it ranges between $1.6$ and $130$.

In addition to the exclusion of the small Gaussian PSF, the difference in the
noise addition and the difference in the ranges of the SNR, in this data set
we also add a few ‘blank’ images (30,000 in all –27,000 for training & 3000
for testing) to the truth data set in contrast to the HW23 data set. Addition
of these blanks will help us determine how the network would perform in
conditions of very low SNR, or when spuriosities in noisy images might appear
object-like even though there is no actual object. All-in-all our data set
consists of $117000$ objects in the training set and $13000$ objects in the
test set including ‘blanks’ and ‘non-blanks’.

##  3 Approach

Convolutional auto-encoders (CAE) are an attractive method for analyzing
astronomical images because they are inherently translationally invariant. One
can train them to perform a deconvolution for a fixed PSF, where the shape of
that PSF gets burned into weights that do the deconvolution. But which
approach should one take if the PSF shape is free and is one of the inputs to
the problem?

Simply feeding the PSF image in addition to the actual noisy input image will
not do the job in an ordinary convolutional network, since they respect
locality: features around a certain coordinate $(x,y)$ are combined in some
complicated non-linear, but still approximately local manner to produce the
output image, while the relation between the PSF and the noisy image is
convolutional:

$I_{i}=I_{gt,i}\circledast PSF_{i}+N_{i},$ (1)

where $I_{i}$ denotes the observed image, $I_{gt}$ the ground truth and
$N_{i}$ the noise image. However, we can let operations that appear in
classical image analysis guide our intuition. These are,

  1. 1.

Convolution with PSF, which we denote as $\times$PSF. Convolution of the noisy
image with PSF acts as a matched filter for object detection and is optimal
for point sources and approximate for small sources. However, it makes any
shape distortion due to assymetric PSF even worse.

  2. 2.

Convolution with inverse PSF, which we denote as $\times$iPSF. Convolution
with the inverse PSF (formally defined as a function which, when convolved
with PSF produces a delta function response) is the naive solution to the
deconvolution problem. When not regularized, it is expected to amplify noise,
but since the standard regularization is to convolve again with a regularizing
kernel (typically a Gaussian), we anticipate that the network would learn this
process by itself.

  3. 3.

Convolution with the inverse PSF convolved with a 1 arcsec Gaussian PSF, which
we denote as $\times$iPSF1. This is a regularized version of the above, which
can make it better, but also looses information. The size of the regularizing
kernel was chosen after several trial and error scenarios, from which it was
concluded that $1$ is the standard deviation at which the network performs the
best. We have confirmed that our results are stable across reasonable range of
regularizing kernel sizes.

For each actual input galaxy we produced these 3 types of manipulated data.
Each image of each of these types has a size of $32\times 32\times 3$ (trimmed
down from the $35\times 35\times 3$ used in HW23 and preprocessed in the same
manner with maximum value scaling). While in HW23 the maximum values across
each object for the noisy and truth data sets was used for scaling, in our
case the maximum values across each object of the $\times$PSF, $\times$iPSF,
$\times$iPSF1and the truth data is used. In this manner, the minimum and
maximum values of all four data sets fall between 0 and 1. We experimented
with various combinations of those inputs to identify those that work
optimally and therefore our network input size varied from $32\times 32\times
3$ to $32\times 32\times 9$.

###  3.1 Convolutional Autoencoder

The network used in this work is a simple 2-dimensional convolutional
autoencoder (CAE) with three layers ($\times 2$ combinedly for the encoder and
the decoder parts).

Autoencoders (AE, Rumelhart et al. 1985) are a “self-supervised” machine
learning method that is generally used in image processing for dimensionality
reduction and feature extraction with the encoder and decoder parts being
separate neural networks with similar configurations. The encoder takes the
input data and maps it to a latent space (encoded space) which is generally of
a lower dimensionality than the input data, thus achieving data compression
and dimensionality reduction. The decoder decompresses/recompresses the data,
often with some associated loss which is dealt with by neural networks in the
normal fashion, such as by using backpropagation to update weights between
batches. The bottleneck thus constructed allows for the clearly structured
portion of the data to pass through, which brings about effective denoising
(an important consideration, especially in astronomical images). In the case
of a single layer AE, this framework is analogous with principal component
analysis (PCA, Pearson 1901; Hotelling 1933) provided it uses a linear
activation function. When the AE consists of multiple layers, thus making it
deep and non-linear, the loss during the decoder phase becomes more complex.
This is the advantage AE has over PCA, it’s ability to learn non-linear
patterns in the input data with lower dimensions and less data loss. As
mentioned earlier, in this work we have implemented a CAE, which is simply an
AE with convolutional layers and a convolutional bottleneck, and in our case,
is trained end-to-end111There are also AEs that are trained layer by layer,
but these are a different type of AEs called “stacked” AEs.. More details
about the specific network architecture is given in the following section.

###  3.2 Network architecture

Our implementation of the CAE is based on the Keras library in Python
222https://keras.io/ within Tensorflow framework
333https://www.tensorflow.org/, utilising the available APIs. The encoder
layer consists of three convolutional layers and three pooling layers with the
convolutional layers activated with a ‘Leaky ReLU’ (Maas et al. 2013, leaky
rectified linear unit) function with $\alpha=0.1$. In normal ReLU (Glorot et
al., 2011), the negative part of the function is set to 0 which renders the
unit inactive. This causes what is called ‘dying ReLU’ problem that can
sometimes lead to overfitting. In Leaky ReLU this is dealt with by applying a
non-zero slope, indicated by $\alpha$. All $6$ layers have padding set to
‘same’ and the pooling utility used for spatial downsampling is ‘Maxpooling’.
The kernels used are $5\times 5$ in the convolutional layers and $2\times 2$
in the pooling layers. The convolutional layers have respective filter sizes
of 64, 128 and 256. The decoder layers are symmetric to the encoder layers
except that a transpose convolution is applied and that instead of the pooling
layers there is an ‘UpSampling’ layer per transpose convolution layer for
spatial upsampling, which also utilises a $2\times 2$ kernel. The final output
layer has 3 units with a $3\times 3$ kernel, ‘same’ padding and activated with
a ‘softplus’ function (Nair & Hinton, 2010). Figure 1 shows a schematic
representation of the network.

Figure 1: Figure showing the schematic representation of the CAE used in this
work. The input image is a compound data cube of $\times$PSF, $\times$iPSFand
$\times$iPSF1with dimensions $32\times 32\times 9$ per image. The output image
dimensions are $32\times 32\times 3$. The numbers shown indicate the change in
dimensions of the image tensor as it propagates through the network. The
different components of the network are indicated in the legend. The network
diagram was generated using Bäuerle et al. (2021).

The regularization aspect of the network is achieved by using
‘BatchNormalization’ layers in both the encoder and decoder portions. There is
much discussion in literature as to whether dropout is necessary for CNNs or
whether batch normalization is useful for autoencoders. In this case, based on
trial and error we have decided to simply apply ‘BatchNormalization’ layers
after each pooling/upsampling layer as might be the case in the final network.
The model thus formulated is compiled with the Adam optimizer (Kingma & Ba,
2014) and the loss function BinaryCrossentropy
444https://www.tensorflow.org/api_docs/python/tf/keras/losses/BinaryCrossentropy,
both from the Keras API. The initial learning rate is $0.0001$ and it decays
exponentially at the rate of $0.6$ per $100000$ steps. The model is trained
for $350$ epochs with a batch size of $10$. We intend for this method to be a
more specific methodological sequel to the work in HW23 in order to provide
unbiased shear measurements that do not correlate with the PSF but do so with
the truth image. This simple, ‘quick and dirty’ approach does indeed provide
results comparable to that of HW23 for specific cases as we will illustrate in
the following sections. It is also cheaper in terms of computation time and
resources, taking $\sim 3$ hours to train on an NVIDIA GeForce GTX 1650 GPU
with 4GB memory and compiled with cuDNN555https://developer.nvidia.com/cudnn
(version 11.4) .

###  3.3 (Re-)Normalization

A perennial problem in the application of neural networks to astronomical
image analysis is the image dynamic range. Dynamic range in astronomical
images can vary by many orders of magnitude, while the neural networks require
inputs to be of limited dynamic range. The non-linear transformations inside
the neural network are typically non-linear over the range $0-1$ and if given
an input range that is very different, the system will fundamentally change
its properties. The usual procedure is therefore to apply certain
normalization factors to the inputs, outputs as well as the truth images. The
main issue is that one is “not allowed” to derive these factors from truth
images, since for real-life problems we do not have truth images. Instead we
argue that the normalization factor should be recalculated with an afterburner
maximum likelihood for the amplitude parameter only.

In our case the procedure is as follows:

  1. 1.

Train the network with input and output images normalized in any sensible
manner, e.g. using min-max rescaling to unity interval.

  2. 2.

When applying to the data, use the same rescaling on the input image as in the
training. The output image is now deconvolved, but normalized in an arbitrary
manner.

  3. 3.

Perform a maximum-likelihood fit to the output image normalization. In our
case, this involves reconvolving the output image with the PSF and then
calculating the renormalization factor as

$r=\frac{(\sum_{{\rm pixels}\ i)}{n_{i}p_{i}}}{\sum_{({\rm pixels}\
i)}{p_{i}^{2}}},$ (2)

where $p_{i}$ is the output image convolved with PSF and $n$ is the noisy
image. This formula can be derived by minimizing the square difference between
the $p$ and $n$.

This approach leverages the neural network to do the actual heavy lifting in
terms of determining the morphology of the output image, while leveraging the
exact solution for the image normalization part. A similar approach can also
be used in various deblending approaches employing neural networks, where
perhaps multiple component amplitudes can be computed.

Upon implementing this procedure, we have found that the measured fluxes have
small negative biases. This is due to the fact that any imperfection in the
recovered shape will result in the object amplitude being biased low. To see
this, imagine if the image to be fit in amplitude has a source at a completely
wrong position. That source would not be able to model the flux of the actual
source and therefore its flux would be low (and presumably consistent with the
upper limit for the flux of a putative source at the wrong position). In our
case, the situation is less dramatic, but sources are still biased low at a
few percent. In order to fix this, we model the true total flux $M_{00}$ as a
quadratic function in the recovered flux. The motivation for the quadratic
function is that we expect morphology to be better recovered with higher
fluxes and therefore less biased.

##  4 Experiments & Results

###  4.1 Deconvolution performance

---
Low SNR

---

Medium SNR

---

High SNR

Figure 2: Figure showing some typical results from the
$\times$PSF-$\times$iPSF-$\times$iPSF1 network. The columns correspond to 1)
noisy image, 2) noisy image convolved with the image PSF ($\times$PSF), 3)
noisy image convolved with the inverse of the image PSF ($\times$iPSF) 4)noisy
image convolved with the image PSF to which a Gaussian filter of standard
deviation $1$ is applied ($\times$iPSF1), 5) truth image and 6) the prediction
from the network. The rows are arranged according to the SNR as indicated.

We start by conducting experiments with different combinations of the
aforementioned $\times$PSF, $\times$iPSFand $\times$iPSF1 data configurations.
Some typical results using the data cube
$\times$PSF-$\times$iPSF-$\times$iPSF1 are represented in Figure 2. As
indicated in the postage stamps, the columns are noisy data, $\times$PSF,
$\times$iPSF, $\times$iPSF1, truth and prediction respectively, arranged by
the image SNR. We would like to note that these are only the results for the
non-blank objects, obtained from a network that was trained with blank images
as well.

We see that the network is capable of successfully deconvolving fine
structures well below the size of the PSF. The $\times$iPSF1 image is close to
a formally regularized solution (e.g. Wiener-filter like). For low SNR objects
in Fig 2 it under-regularizes, as indicated by still significant noise
features, while for the high SNR it over-regularizes. In any case, the
prediction can correctly deconvolve fine structures in the truth image even in
the low-noise regime. These images show that as the basic level, the network
performs as one would expect. We next turn to more quantitative assessments of
its performance.

|  $\times$PSF-$\times$iPSF-$\times$iPSF1 |  $\times$PSF-$\times$iPSF |  $\times$iPSF-$\times$iPSF1 |  $\times$PSF-$\times$iPSF1 |  $\times$PSF |  $\times$iPSF |  $\times$iPSF1
---|---|---|---|---|---|---|---
| Mean | Median | Mean | Median | Mean | Median | Mean | Median | Mean | Median | Mean | Median | Mean | Median
PSNR | 75.177 | 74.744 | 73.477 | 72.873 | 74.056 | 73.611 | 75.378 | 74.987 | 74.034 | 73.709 | 73.331 | 72.842 | 75.290 | 74.943
1-SSIM | 0.0489 | 0.0177 | 0.0557 | 0.0253 | 0.0507 | 0.0209 | 0.0404 | 0.0165 | 0.0628 | 0.0259 | 0.0708 | 0.0281 | 0.0418 | 0.0173

Table 1: Table showing the differences between the mean and median in the PSNR
(dB) and SSIM values for the truth and predicted data sets for the different
data cubes. For both PSNR & SSIM, the higher the value, the better the
congruence between the truth and the prediction.

In Table 1, we compare the mean and median values of certain standard metrics
that are typically used in image analysis for all the data combinations. Peak
signal-to-noise ratio (PSNR) is defined as the ratio of maximum possible
signal strength to the distorting noise level, expressed in logarithmic
decibel scale. It performs a pixel-to-pixel comparison between images, and
higher values of this quantity are equated to better quality, and lower values
to greater numerical dissimilarities between images (Horé & Ziou 2010). SSIM
on the other hand, defines the perceived similarity between images through the
correlation of image pixels, thus providing information about contrast,
luminance and structure. It maps the structural similarities between two
objects, and is often taken as an acceptable proxy to human visual perception,
and is as such, considered to be more sensitive to image degradation as a
result of compression, noise etc. than PSNR. Here we have used the
implementations in scikit-image (Van der Walt et al. 2014) to calculate both
PSNR and SSIM. In general, the higher the PSNR and SSIM values are, the
better. Here, for SSIM we have opted to show how dissimilar the truth and
prediction images are (1-SSIM since the maximum possible value is 1) rather
than the values themselves, so the lower the value the better the similarity
between the truth and prediction. For both PSNR and SSIM, the best mean and
median values is returned by $\times$PSF-$\times$iPSF1, closely followed by
the full $\times$PSF-$\times$iPSF-$\times$iPSF1 while $\times$iPSF returns the
worst mean and median values.

|  $\times$PSF-$\times$iPSF-$\times$iPSF1 |  $\times$PSF-$\times$iPSF |  $\times$iPSF-$\times$iPSF1 |  $\times$PSF-$\times$iPSF1 |  $\times$PSF |  $\times$iPSF |  $\times$iPSF1
---|---|---|---|---|---|---|---
| Mean | RMS | Mean | RMS | Mean | RMS | Mean | RMS | Mean | RMS | Mean | RMS | Mean | RMS
$M_{00}$ | 0.4313 | 17.6943 | 0.8253 | 24.4059 | 0.6633 | 24.5903 | 0.9567 | 17.9097 | -0.5118 | 26.4617 | 0.1937 | 26.2412 | 0.9077 | 18.5814
$\left<x\right>$ | 0.0044 | 0.6096 | 0.0293 | 0.6314 | 0.0201 | 0.6543 | 0.0392 | 0.5610 | 0.0237 | 0.5391 | 0.0192 | 0.8376 | 0.0337 | 0.6029
$\left<y\right>$ | 0.0070 | 0.6034 | -0.0166 | 0.6204 | 0.0033 | 0.6282 | -0.0070 | 0.5436 | -0.0229 | 0.5486 | 0.0365 | 0.863 | 0.0144 | 0.5575
$e_{1}$ | -0.0005 | 0.1443 | -0.0050 | 0.1503 | -0.0082 | 0.1424 | -0.0030 | 0.1376 | 0.0108 | 0.1541 | -0.0043 | 0.1585 | 0.0010 | 0.1350
$e_{2}$ | 0.0042 | 0.1586 | -0.0147 | 0.1719 | 0.0040 | 0.1618 | -0.0011 | 0.1561 | -0.0117 | 0.1717 | -0.1120 | 0.1862 | -0.0114 | 0.1605
$\|\boldsymbol{e}\|$ | 0.0670 | 0.1425 | 0.0581 | 0.1465 | 0.0486 | 0.1356 | 0.0444 | 0.1308 | 0.0908 | 0.1698 | 0.0729 | 0.1581 | 0.0425 | 0.1281
$s$ | -4.3078 | 11.9458 | -2.8706 | 12.2596 | -2.4115 | 12.2552 | -0.3480 | 9.2682 | -8.7653 | 16.4317 | -6.7513 | 18.1492 | -0.8512 | 9.1547

Table 2: Table showing the mean and RMS values of the differences in moments &
other astronomically significant quantities between the true and the predicted
images for the different data combinations.

Table 2 shows a comparison of quantities that are physically significant in
astronomical image processing in terms of image moments (see Equation 10 in
HW23). We have tried to quantify our results in terms of the four major image
aspects, i.e., (1) position, (2) flux, (3) shape and (4) size. Here, the
central moment $M_{00}$ is a measure of the total flux, $<x>$ & $<y>$
astrometric co-ordinates ($M_{10}$ and $M_{01}$ moments), the second central
moments $x^{2}$, & $y^{2}$ ($\mu_{20}$ and $\mu_{02}$) representing size,
$s=x^{2}+y^{2}$, and $e_{1}$ & $e_{2}$ the $+$ and $-$ components of
ellipticity, and $\|\boldsymbol{e}\|=e_{1}^{2}+e_{2}^{2}$. Since these are the
differences between the mean and RMS values of the truth and prediction, the
objective is that the absolute value of the difference be as small as
possible.

Figure 3: 2-D histogram showing the $M_{00}$ moments of the raw truth data
(X-axis) and the re-scaled and flux corrected predictions (Y-axis). We use
this to illustrate the accuracy of our re-scaling method and flux corrections.

We apply corrections to the recovered $M_{00}$ values as described in the
Section 3.3. The resulting function for the
$\times$PSF-$\times$iPSF-$\times$iPSF1 has linear and quadratic coefficients
equal to $0.984$ and $2.84\times 10^{-6}$ and a constant term $-1.6$. The
error percentage or bias between the predictions and this corrected data are
$7.25\%$, $2.15\%$ and $1.5\%$ respectively for the low, medium and high SNR
bins and those between the truth data and the corrected predictions are
$1.46\%$, $0.22\%$ and $0.066\%$, well within the noise scatter and consistent
with no bias.

In Figure 3, we illustrate these results graphically, We plot a 2-d histogram
of the $M_{00}$ values for the raw truth data and the prediction values
rescaled using the factors calculated in Equation 2 and then flux corrected
for the $\times$PSF-$\times$iPSF-$\times$iPSF1 data cube. As can be seen, the
$M_{00}$s are within the same range and to a great extent align with each
other. We single out the $M_{00}$ moment here to illustrate this effect since
all the other quantities that we consider are flux scale independent. Even
though only one of the metrics we consider is affected by it, the impact of
re-scaling cannot be ignored because $M_{00}$ is representative of the total
flux and is therefore quite important in terms of image recovery. We note that
all fluxes are biased low. That is because any shape mismatch between the
output image and the truth will result in normalization being low.

The best values for both the mean $M_{00}$ (denoting bias) is returned by the
$\times$iPSF1 data cube, and the RMS (scatter) by
$\times$PSF-$\times$iPSF-$\times$iPSF1, the worst mean by
$\times$PSF-$\times$iPSF1, and the worst RMS by $\times$PSF with the other
data configurations performing moderately. The best mean value for
$\left<x\right>$ is returned by $\times$PSF-$\times$iPSF-$\times$iPSF1 and
that for $\left<y\right>$ by $\times$iPSF-$\times$iPSF1. The latter is quite
interesting because convolving noisy image with the image PSF ($\times$PSF) is
a matched filter for point sources and therefore often used for object
detection in astronomical images, but here the object detection aspect seems
to be achieved better (albeit marginally) without the $\times$PSF component.
It is also worth noting that for mean $\left<y\right>$,
$\times$PSF-$\times$iPSF-$\times$iPSF1 and $\times$PSF-$\times$iPSF1 have
remarkably similar values while for $\left<x\right>$ the values of
$\times$PSF-$\times$iPSF-$\times$iPSF1 is an order of magnitude better than
all others. For the RMS value of the astrometric co-ordinates however,
$\times$PSF returns the best (lowest) values. All other data cubes exhibit
similar performances except for $\times$iPSF which performs the worst.

The best mean value for the $+$ component of ellipticity, $e_{1}$, is returned
by $\times$PSF-$\times$iPSF-$\times$iPSF1 and the worst by $\times$PSF;
likewise the best RMS value is returned by $\times$iPSF1, closely followed by
$\times$PSF-$\times$iPSF1. For the $-$ component of ellipticity, $e_{2}$, the
best mean value is returned by $\times$iPSF-$\times$iPSF1 &
$\times$PSF-$\times$iPSF-$\times$iPSF1 (the difference between them only
0.0002) and the worst by xiPSF. The latter is a similar result to that of the
$\times$PSF being inept at determining object position, because convolution
with an inverse PSF ($\times$iPSF) is generally implemented for shape recovery
in images. The best RMS value for $e_{2}$ was also returned by
$\times$iPSF-$\times$iPSF1 and the worst by both $\times$PSF-$\times$iPSF &
$\times$PSF (also 0.0002 difference). For the aggregated ellipticity,
$|\boldsymbol{e}\|$, the best mean and RMS values are returned by
$\times$iPSF1 and the worst of both by $\times$PSF. For the quantity $s$,
denoting the sizes, the best mean value is returned by
$\times$PSF-$\times$iPSF, the best RMS by $\times$iPSF1, the worst mean by
$\times$PSF and the worst RMS by $\times$iPSF.

In Figure 4, we explore the output ellipticity correlations with PSF & truth
image using Equations (29) and (30) in HW23 (Eqs. 29 & 30 hereafter). The left
panel denotes the PSF correlation of the predicted and the truth data sets for
the different experiments that we have conducted. For an ideal deconvolution,
this correlation should be zero.

This is an important quantity since in applications like shear estimation,
presence of residual PSF dependence could lead to ‘additive bias’ where shapes
are contaminated by PSF (see Mandelbaum 2018 and references therein). An ideal
network should remove all PSF dependence in the predicted image, leaving the
quantity defined in Eq. 29, $\xi_{PSF}=0$. In order to verify this, we also
examined the recovered shapes of objects in relation with the truth image as
defined in Eq. 30. For an ideal case this quantity, $\xi_{true}=1$ and this is
represented in the right panel of the figure. The different data
configurations used as input to the network are denoted in the legend. As is
evident from the left panel of the figure, it seems that just $\times$PSF and
$\times$PSF-$\times$iPSF both do not seem to provide sufficient and
information for the network to remove the residual PSF dependence. All the
other combinations seem to be comparable in this respect. We expect the
$\times$PSF case to be the worst since the information is simply not present
and we indeed find this to be the case. The & $\times$PSF-$\times$iPSFalso
does quite bad, likely by being dominated by the information present in the
$\times$PSF. All the data combinations have similar performance with regard to
ellipticity recovery as can be seen from the right panel of Figure 4 with
marginally better performances by $\times$PSF-$\times$iPSF1 in the low SNR
bin, and $\times$iPSF1 in the medium and high SNR bins. For convenience of
visualization, the points are slightly offset on the X-axis.

Figure 4: Plot showing the residual PSF dependence (left) and correlation with
the truth data (right), separated by SNR. The different data combinations are
indicated in the legend below, and they are slightly offset in the X-axis for
better visualization.

###  4.2 Object detection efficiency

Figure 5 shows the confusion matrices for a few selected networks, starting
with $\times$PSF (our worst performing network) to the complete
$\times$PSF-$\times$iPSF-$\times$iPSF1 network. Confusion matrices are a
standard representation method for the results of machine learning algorithms,
especially classification problems. In a typical $2\times 2$ confusion matrix
used for binary classification problems, the values represented are the true
positive rate (TPR; predictions and truth correspond), false positive rate
(FPR; predictions classify as positive while the truth disagrees) from left to
right in the first row; the false negative rate (FNR; predictions classify as
negative while the truth disagrees) and the true negative rate (TNR;
predictions and truth agree that it’s negative) from left to right in the
second row. The left diagonal represents the correct classifications. In our
case, we can divide our data set into four classes – blanks, low SNR, medium
SNR and high SNR. In reconstruction problems such as ours, a threshold needs
to be defined to classify the predictions into detections or non-detections.
We defined these thresholds upon visual inspection of the histogram
representations of the SNRs of the flux corrected predictions.

(a) $\times$PSF (b) $\times$iPSF1

(c) $\times$PSF-$\times$iPSF (d) $\times$PSF-$\times$iPSF-$\times$iPSF1

Figure 5: Figure showing the confusion matrices for a few selected networks.
The true labels are on the Y-axis and the predicted labels on the X-axis.

Part (a) of Figure 5 represents the confusion matrix for the network trained
solely using $\times$PSF. For the blanks bin, this network seems to return the
best results, which makes sense intuitively, since noisy images are often
convolved with their respective PSFs for object detection, and when there is
no object, this must make the task easier. However, this reasoning breaks down
at the threshold between the blank and low bins for this network and between
low and medium bins for this network. The blank vs low confusion seems to be
prevalent in the case of all the different networks with varying degrees, it
seems especially worse in the case of $\times$PSF. The network trained with
$\times$PSF-$\times$iPSF(part (c) seems to perform the best in this respect
with an accuracy of $77\%$. The differentiation between the low and the medium
bins seems to improve with addition of more/slightly diverse information, as
can be seen with the other networks – $\times$iPSF1 in particular seems to
deal with this better. The best performance in the medium bin is returned by
$\times$iPSF1(part (b)), with an accuracy of $\sim 99\%$. The other networks
also return comparable performance TPRs in this bin, indicating an expected
trend of the accuracy improving and stabilising as the SNR increases. Part (c)
of the figure shows the results from the network trained on
$\times$PSF-$\times$iPSF-$\times$iPSF1, and it returns the a marginally high
accuracy compared to the other networks in the high SNR bin, $\sim 99\%$.

From both Figure 5 and from Tables 1–2 we can comfortably conclude that the
inclusion of $\times$iPSF1, while unconventional, is definitely helpful in
both detection and shape recovery, especially in the low SNR regime. It is
also similarly evident that only $\times$PSF or only $\times$iPSF is simply
not sufficient information for this purpose. As can be seen from Figure 2,
even visually it appears that $\times$iPSF1 encompasses salient features of
both $\times$PSF and $\times$iPSF, which could be the reason why it is more
suited to the purpose at hand. It is also possible that this is data dependent
– while we have tried to make our data set as varied as possible, this method
might perform differently in a distinct data set, requiring certain
adjustments, for e.g. the Gaussian regularising filter might have a different
FWHM. The size of the data cube is not linearly related to the running time
required, either way the network takes between 2 and 3 hours to complete,
which provides a solid scope for experimentation with the data manipulation
parameters.

##  5 Summary & arbitrary conclusions

We present in this work, a simpler approach to remove residual PSF dependence
and ellipticity recovery as a follow-up to the work done in HW23. The network
implemented in this work is a convolutional autoencoder with 3 encoder and
correspondingly 3 decoder layers with respective filter sizes of 64, 128 and
256. $5\times 5$ kernels are used in the convolutional layers (both encoder
and decoder) while $2\times 2$ kernels are used in the pooling layers in the
encoder and upsampling layers in the decoder. All convolutional layers are
activated with the LeakyReLu function and the final output layer with a
softplus function. The model is compiled with the Adam optimiser and the loss
function used is BinaryCrossentropy. We trained this network with several data
configurations, namely, $\times$PSF, $\times$iPSF, $\times$iPSF1 and
combinations thereof, results of which are detailed in the preceding sections.

We quantify our results using different metrics, which are displayed in
Figures 2 \- 5 and Tables 1 \- 2. In Figure 2, we show some typical results
that is output by the network that is trained with a
$\times$PSF-$\times$iPSF-$\times$iPSF1 data cube/combination as a proof of
concept. We analyse this further with specific correlations between the
predicted data and the PSF and between the predicted data and the true
ellipticities in Figure 4 in the three SNR bins. In Tables 1 and 2, we compare
significant quantities, both from image analysis and astronomical perspectives
to quantify the similarities between the predictions and the truth data. In
Table 1, the quantities compared are PSNR and SSIM which are popular metrics
used to evaluate image reconstruction. While PSNR and SSIM provide some
insight into image recovery from an image analysis standpoint, we are more
interested in the recovery of astronomically significant features. To this
end, in Table 2, we compare the differences between moment proxies for total
fluxes, first order and second order moments averaged over the three bands for
the different data combinations. Then we analysed our results using confusion
matrices, a very popular form of representation generally used in
classification problems by defining thresholds based on SNR, as shown in
Figure 5.

We conclude that as is general convention, just convolving the noisy image
with the PSF ($\times$PSF) or with the inverse PSF ($\times$iPSF) by
themselves are not informative enough for the autoencoder network to
efficiently locate the object and remove the PSF dependence. $\times$iPSF1 on
the other hand, provides very good results, both by itself and when used in
combination with $\times$PSF or $\times$iPSF1 (or both).

It is hard to declare a clear winner. $\times$iPSF1 is the most economical and
is close to the best, although on certain tests, augmenting it with the other
two combinations improves results. The simplest interpretation is that
$\times$iPSF1 provides a blurry image that is standardized across various
observed PSFs – the role of the network is then to simply sharpen and denoise
it based on how galaxies in the training set look like. As in any denoising /
deblurring approach, the missing information is filled-in based on priors and
therefore one should be cautious not to over-interpret galaxy morphologies
derived this way.

The main advantages that we see for this method that we have put forth are,
(1) simple neural network (2) fast processing, (3) minimal data manipulation
and (4) minimal reliance on data features that might not be readily available
(e.g. normalization factors). We would also like to emphasize that this method
is not a ‘be all and end all’ product in removing PSF dependence and
ellipticity recovery by any means. This is merely a simpler way to obtain fast
results in this respect using marginal data manipulation. It is possible that
its performance will vary in a different data set. While we haven’t tested it
on other data sets, we feel that it is not suitable for transfer learning
purposes. Although that requirement might be rendered moot by the fact that
this is a very fast method even on a data set that contains some 100,000
objects, so training on a similar sized data set would also be a fast process.
Our only claim is that for similar data that has both noise and contains one
centred object, this is an effective, albeit crude, specific method to obtain
denoised data that has excellent correlation to the truth image. There are
certain approaches that might be applied to refine this method further. For
example, it might be interesting to replace the CAE with a VAE – which would
add both some probabilistic and Bayesian aspects to the denoising process.
Moreover, VAEs also provide more tuneable parameters that give more control
over the latent representation of the input data.

## References

  * Abbott et al. (2018) Abbott T. M. C., et al., 2018, Phys. Rev. D, 98, 043526
  * Arcelin et al. (2021) Arcelin B., Doux C., Aubourg E., Roucelle C., LSST Dark Energy Science Collaboration 2021, MNRAS, 500, 531
  * Bertin & Arnouts (1996) Bertin E., Arnouts S., 1996, Astronomy and astrophysics supplement series, 117, 393
  * Boucaud et al. (2020) Boucaud A., et al., 2020, Monthly Notices of the Royal Astronomical Society, 491, 2481
  * Bäuerle et al. (2021) Bäuerle A., van Onzenoodt C., Ropinski T., 2021, IEEE Transactions on Visualization and Computer Graphics, 27, 2980
  * Gardner et al. (2006) Gardner J. P., et al., 2006, Space Sci. Rev., 123, 485
  * Glorot et al. (2011) Glorot X., Bordes A., Bengio Y., 2011, in Proceedings of the fourteenth international conference on artificial intelligence and statistics. pp 315–323
  * Green et al. (2012) Green J., et al., 2012, arXiv e-prints, p. arXiv:1208.4012
  * Hikage et al. (2019) Hikage C., et al., 2019, PASJ, 71, 43
  * Horé & Ziou (2010) Horé A., Ziou D., 2010, in 2010 20th International Conference on Pattern Recognition. pp 2366–2369, doi:10.1109/ICPR.2010.579
  * Hotelling (1933) Hotelling H., 1933, Journal of educational psychology, 24, 417
  * Kingma & Ba (2014) Kingma D. P., Ba J., 2014, arXiv preprint arXiv:1412.6980
  * Kingma & Welling (2013) Kingma D. P., Welling M., 2013, arXiv preprint arXiv:1312.6114
  * Koekemoer et al. (2007) Koekemoer A. M., et al., 2007, The Astrophysical Journal Supplement Series, 172, 196
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration et al., 2009, arXiv e-prints, p. arXiv:0912.0201
  * Lanusse et al. (2021) Lanusse F., Mandelbaum R., Ravanbakhsh S., Li C.-L., Freeman P., Póczos B., 2021, Monthly Notices of the Royal Astronomical Society, 504, 5543
  * Maas et al. (2013) Maas A. L., Hannun A. Y., Ng A. Y., et al., 2013, in Proc. icml. p. 3
  * Mandelbaum (2018) Mandelbaum R., 2018, ARA&A, 56, 393
  * Nair & Hinton (2010) Nair V., Hinton G. E., 2010, in Proceedings of the 27th international conference on machine learning (ICML-10). pp 807–814
  * Pearson (1901) Pearson K., 1901, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 6, 559
  * Refregier et al. (2010) Refregier A., Amara A., Kitching T. D., Rassat A., Scaramella R., Weller J., 2010, arXiv e-prints, p. arXiv:1001.0061
  * Reiman & Göhre (2019) Reiman D. M., Göhre B. E., 2019, Monthly Notices of the Royal Astronomical Society, 485
  * Ronneberger et al. (2015) Ronneberger O., Fischer P., Brox T., 2015, arXiv e-prints, p. arXiv:1505.04597
  * Rowe et al. (2015) Rowe B., et al., 2015, Astronomy and Computing, 10, 121
  * Rumelhart et al. (1985) Rumelhart D. E., Hinton G. E., Williams R. J., et al., 1985, Learning internal representations by error propagation
  * Scoville et al. (2007) Scoville N., et al., 2007, The Astrophysical Journal Supplement Series, 172, 38
  * Van der Walt et al. (2014) Van der Walt S., Schönberger J. L., Nunez-Iglesias J., Boulogne F., Warner J. D., Yager N., Gouillart E., Yu T., 2014, PeerJ, 2, e453
  * Wang et al. (2022) Wang H., Sreejith S., Slosar A. c. v., Lin Y., Yoo S., 2022, Phys. Rev. D, 106, 063023
  * Wang et al. (2023) Wang H., Sreejith S., Lin Y., Ramachandra N., Solsar A., Yoo S., 2023, The Open Journal of Astrophysics, 6, 30


#  From asymptotic freedom to $\theta$ vacua:
Qubit embeddings of the O(3) nonlinear $\sigma$ model

Stephan Caspar  0000-0002-3658-9158<EMAIL_ADDRESS>Hersh Singh
0000-0002-2002-6959<EMAIL_ADDRESS>InQubator for Quantum Simulation (IQuS),
Department of Physics, University of Washington, Seattle, Washington
98195-1550, USA

###### Abstract

Conventional lattice formulations of $\theta$ vacua in the $1+1$-dimensional
$\text{O}(3)$ nonlinear sigma model suffer from a sign problem. Here, we
construct the first sign-problem-free regularization for _arbitrary_ $\theta$.
Using efficient lattice Monte Carlo algorithms, we demonstrate how a
Hamiltonian model of spin-$\tfrac{1}{2}$ degrees of freedom on a 2-dimensional
spatial lattice reproduces both the infrared sector for arbitrary $\theta$, as
well as the ultraviolet physics of asymptotic freedom. Furthermore, as a model
of qubits on a two-dimensional square lattice with only nearest-neighbor
interactions, it is naturally suited for studying the physics of $\theta$
vacua and asymptotic freedom on near-term quantum devices. Our construction
generalizes to $\theta$ vacua in all $\text{CP}(N-1)$ models, solving a long
standing sign problem.

††preprint: IQuS@UW-21-025, INT-PUB-22-012

##  I Introduction

The strong interactions of the standard model described by quantum
chromodynamics (qcd) pose a challenging problem for classical computation.
While nonperturbative lattice Monte Carlo (mc) methods are a powerful tool for
studying static properties of strongly coupled quantum field theories (qfts)
like qcd [1, 2, 3, 4], questions involving real time dynamics, finite density
or nontrivial $\theta$ vacua are still out of reach for lattice mc methods due
to severe sign problems [5, 6].

Emerging quantum platforms provide an exciting possibility for investigating
qfts in previously inaccessible regimes. They are not directly affected by the
sign problems arising in classical lattice mc methods. However, bosonic
lattice field theories such as qcd have infinite-dimensional local Hilbert
spaces, while hardware degrees of freedom (dof) are usually finite-
dimensional, mostly qubits. A significant effort is underway to explore
different embeddings of qfts in qubits, with a multitude of ideas emerging
from bosonic field theory [7, 8, 9, 10], nonlinear sigma models (nl$\sigma$ms)
[11, 12, 13, 14, 15, 16, 17, 18, 19, 20] and gauge theories [21, 22, 23, 24,
25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37].

The $1+1$-dimensional $\text{O}(3)$ nl$\sigma$m has a long history as a
prototype for qcd, due to similarities such as asymptotic freedom, dynamical
transmutation and the generation of a nonperturbative mass gap, as well as a
topological $\theta$ term. The $\text{O}(3)$ nl$\sigma$m with a $\theta$-term
is formally defined by the continuum action

$\displaystyle S_{\theta}[\vec{\phi}]=\frac{1}{g^{2}}\int
d^{2}x(\partial_{\mu}\vec{\phi})^{2}+i\theta Q[\vec{\phi}],$ (1)

where $\vec{\phi}\in{\mathbb{R}}^{3}$ with $|\vec{\phi}|^{2}=1$, and

$\displaystyle Q[\vec{\phi}]=\frac{1}{8\pi}\int d^{2}x\ \varepsilon_{\mu\nu}\
\vec{\phi}\cdot(\partial^{\mu}\vec{\phi})\times(\partial^{\nu}\vec{\phi})$ (2)

is the integer topological charge, making the theory $2\pi$-periodic in
$\theta$. Both $\theta=0,\pi$ points are well-understood, analytically as well
as on the lattice. Exact S-matrices have been conjectured for both $\theta=0$
and $\theta=\pi$ [38, 39, 40, 41] and their integrability has been confirmed
using non-perturbative lattice mc methods [42, 43, 44].

Figure 1: rg flow diagram of $\text{O}(3)$ nl$\sigma$ms $S_{\theta}$ defined
in Eq. 1. $S_{\theta}$ is a family of asymptotically-free qfts which all flow
into the trivial ir fixed point, except at $\theta=\pi$ where it reaches the
$\text{SU}(2)_{1}$ wzw fixed point. At small $|\theta-\pi|$, the rg flow of
$S_{\theta}$ passes arbitrarily close to the wzw fixed point, on its way to
the trivial fixed point.

However, general, non-integrable $\theta$ remain challenging. As a topological
effect, it cannot be studied directly in perturbation theory about the free
ultraviolet (uv) fixed point, although some analytic progress has been made by
perturbing about the $\theta=\pi$ integrable point [45]. Nonperturbatively,
the inclusion of a $\theta$ term causes a sign problem when discretizing the
action in Eq. 1 on a $2$-dimensional spacetime lattice. Even though improved
actions combined with cluster algorithms have been shown to tame both cutoff
effects and the sign problem to allow a reliable extrapolation from modest
volumes around $\theta\approx 0$ [46, 47], and even fully solve the sign
problem at $\theta=\pi$ [43], so far there are no known lattice mc methods
which allow a fully controlled study of arbitrary $\theta$ vacua.

Motivated by the prospect of quantum simulation to address these challenges,
we develop an embedding of the $\text{O}(3)$ nl$\sigma$m at _arbitrary_
$\theta$ into a $2$-dimensional Heisenberg antiferromagnet, such that a
controlled continuum limit can be taken. Remarkably, not only does this model
allow the systematic study of $\theta$ vacua on quantum hardware, it also
enables the first sign problem free algorithm for classical computations at
arbitrary $\theta$. This extends a similar proposals put forward in Refs. [11,
21, 12, 14] for classical and quantum simulation of $\theta=0,\pi$ theories.

Regularizing qfts using explicitly finite-dimensional local dof is a promising
approach for quantum simulation. Universality lets us understand this
remarkable variety in models with the ability to describe the same continuum
qfts. In his seminal work on rg, Wilson showed how continuum qfts emerge at
second-order critical points of lattice models [48, 49, 50]. In this
framework, the infinite-dimensional continuum fields can arise naturally at
long-distances from finite-dimensional microscopic local dof. While this
approach is natural in the context of quantum computation, universality has
even been leveraged to circumvent sign problems that plague conventional
lattice regularizations. This was shown, for example, with the $\text{O}(3)$
model at finite density [11] and the $\text{CP}(2)$ model at $\theta=\pi$
[51]. Efficient cluster algorithms for $\text{CP}(N-1)$ models have been
demonstrated [13], where a no-go theorem prevents efficient cluster algorithms
using the standard lattice action [52].

Figure 2: Two configurations for the staggered interactions, described in Eq.
4, considered as a regularization of the $1+1$-dimensional $\text{O}(3)$
nl$\sigma$m with a $\theta$ term. For the alternating staggering, all
couplings are antiferromagnetic, while for the columnar case, the transverse
coupling $J^{\prime}$ is ferromagnetic and $J_{\pm}=J(1\pm\gamma)$ is
antiferromagnetic. All interactions are of the Heisenberg
$\vec{S}_{i}\cdot\vec{S}_{j}$ type.

##  II The Qubit Hamiltonian

In this work, we show that the continuum limit of the $1+1$d $\text{O}(3)$
nl$\sigma$m with a $\theta$ term can be obtained from a spin-$\frac{1}{2}$
Heisenberg antiferromagnet on a $2d$ lattice with staggered couplings

$\displaystyle H=\sum_{(x,y)}J_{x,y}\
\vec{S}_{x,y}\cdot\vec{S}_{x+1,y}+J^{\prime}\sum_{(x,y)}\
\vec{S}_{x,y}\cdot\vec{S}_{x,y+1}.$ (3)

where $\vec{S}_{i}$ are the spin operators acting on two-dimensional Hilbert
space at the site $(x,y)$, $J_{x,y}$ are the couplings along the $x$
direction, $J^{\prime}$ is the coupling along the $y$ direction, and the $2d$
lattice has dimensions $L_{X}\times L_{Y}$. We consider the following two
configurations for staggering the couplings:

Alternating: $\displaystyle J^{\prime}>0,\quad J_{x,y}$
$\displaystyle=J(1+(-1)^{x+y}\gamma),$ (4) Columnar: $\displaystyle
J^{\prime}<0,\quad J_{x,y}$ $\displaystyle=J(1+(-1)^{x}\gamma),$

where $J>0$ is always antiferromagnetic, and $\gamma$ is the staggering
parameter, as shown in Fig. 2. In both these cases, the continuum limit of the
$\text{O}(3)$ nl$\sigma$m with a $\theta$ term can be obtained from odd or
even $L_{Y}$, by taking the limit $L_{Y}\to\infty$ at fixed $\gamma L_{Y}$
such that $L_{X}\gg L_{Y}\gg 1$ is maintained.

To demonstrate the continuum limit, we need to recover the physics of the
theory described by Eq. 1 at all scales, from the uv to the ir. For all
$\theta$, the continuum action $S_{\theta}$, defined in Eq. 1, describes an
asymptotically-free theory, controlled in the uv by the fixed point of two
free bosons. The coupling $g$ is a relevant coupling and thus drives the
theory away from the free uv fixed point into a strongly coupled theory in the
ir. While all $S_{\theta}$ theories flow out of the same uv fixed point, non-
perturbative effects lead to different rg trajectories for different $\theta$.
Figure 1 shows a conjectured rg flow diagram for the $\text{O}(3)$ nl$\sigma$m
at arbitrary $0\leq\theta\leq\pi$. For all $\theta\neq\pi$ the theory flows to
the trivial massive fixed point in the ir. However, at $\theta=\pi$, the
theory undergoes a second order phase transition and the low-energy physics
changes completely. The mass-gap vanishes and the ir physics is described by a
nontrivial conformal field theory (cft) called the $\text{SU}(2)_{1}$ wzw
theory [53]. Interestingly, the two ideas of staggering [54, 55] and D-theory
[12] can be combined with the qubit Hamiltonian of Eq. 3 to reproduce the
physics of both ir and uv.

Figure 3: Proposed embedding of the $\text{O}(3)$ nl$\sigma$m with a $\theta$
term into a 2-dimensional array of ultracold atoms. The alternating staggering
described in Eq. 4 and Fig. 2 arises naturally from distance-dependent
antiferromagnetic interactions by deforming a rectangular lattice.

_IR physics of the $\theta$ vacua_. For $\gamma=0$, the Hamiltonian of Eq. 3
reduces to the ordinary Heisenberg antiferromagnet. For a fixed $L_{Y}$ and
$\gamma=0$, this model has been studied in condensed matter literature as spin
ladders, and is known to be described by the $\text{O}(3)$ nl$\sigma$m at low-
energies with $\theta=2\pi SL_{Y}$ [56, 57, 58]. Under this identification,
the translation-by-one symmetry of spins
($\vec{S}_{x,y}\mapsto\vec{S}_{x+1,y}$) on the lattice scale becomes the
charge conjugation symmetry ($\vec{\phi}\mapsto-\vec{\phi}$) in the continuum.
Therefore, a $\theta$-term can be induced in the ir by introducing a staggered
coupling which breaks this symmetry [55, 59, 54]. Ref. [55] showed that for
spin-$S$ ladders with alternating staggering $\theta=2\pi SL_{Y}(1+\gamma
f(L_{Y}))$, where $f(L_{Y})$ is a non-universal function. Therefore the low-
energy physics of the $\theta$ vacua can be studied by varying the staggering
parameter $\gamma$ [60]. However, to obtain the continuum limit of the
$\text{O}(3)$ nl$\sigma$m, we must also obtain the physics of asymptotic
freedom in the uv, which we now turn to.

_Regulating the UV: Asymptotic Freedom._ The continuum limit of $S_{\theta}$,
in Eq. 1, can be obtained from this Hamiltonian model by considering the limit
$L_{Y}\to\infty$ while maintaining $L_{X}\gg L_{Y}\gg 1$. This approach has
been developed under the name D-theory [11, 12, 21], and works as follows. In
the thermodynamic limit $L_{X},L_{Y}\to\infty$, the ground state of the
Heisenberg antiferromagnet has Néel ordering with spontaneously broken global
$SU(2)$-symmetry, with massless Goldstone mode excitations. As $L_{Y}$ is made
finite, the system develops an exponentially large correlation length $\xi\sim
e^{\\#L_{Y}}\gg L_{Y}$. The physics is therefore frozen along the $L_{Y}$
direction, and the system effectively becomes one-dimensional, described by
the $1+1$d $\text{O}(3)$ nl$\sigma$m with an effective coupling $g^{2}\sim
1/L_{Y}$. Since the correlation length diverges exponentially in $L_{Y}$, a
continuum qft can be by defined in the limit of $L_{Y}$ large. Therefore, in
this limit, the spin-$\frac{1}{2}$ Hamiltonian of Eq. 3 is a lattice
regularization of the $\text{O}(3)$ nl$\sigma$m with an arbitrary $\theta$ at
all scales, including asymptotic freedom in the uv.

_Extension to $\text{CP}(N-1)$ models._ All methods in this paper are
straightforward to extend from $\text{O}(3)=\text{CP}(1)$ to the entire family
of $\text{CP}(N-1)$ models, which also allow for a $\theta$ term. Both
$\theta=0,\pi$ have been considered before in the D-theory formulation [14,
15, 51, 13] using a Heisenberg model of $\text{SU}(N)$ spins, where the
$\text{SU}(N)$ representations are chosen such that spontaneous symmetry
breaking (ssb) of the type $\text{SU}(N)\to\text{U}(N-1)$ occurs [61, 62].
This ensures that the continuum $\text{CP}(N-1)$ fields arise as Goldstone
modes as the continuum limit ($L_{Y}\to\infty$) is taken. Since the discussion
of charge conjugation symmetry is identical to that for the $\text{O}(3)$
model, the staggering patterns of Eq. 4 will induce $\theta\neq 0,\pi$ in
these constructions as well.

Figure 4: Step-scaling function of the $\text{O}(3)$ nl$\sigma$m with various
$\theta$. We show step-scaling curves for different values of $\gamma
L_{Y}\sim|\theta-\pi|/\pi$ with odd $L_{Y}$, obtained from alternating (left)
and columnar staggering (right), as defined in Eq. 4. For a fixed $\gamma
L_{Y}$, we show mc results for $L_{Y}=3$ (solid line), $L_{Y}=5$ (dashed
lines) and $L_{Y}=7$ (dotted lines). The dotted black curve is a two-loop
perturbative prediction [44]. The dashed black line is the step-scaling
function for $\theta=0$ obtained in Ref. [44]. The solid black line is an
$O(z^{-10})$ fit to the $\gamma L_{Y}=0$ data, which corresponds to the step-
scaling function of the $\text{O}(3)$ nl$\sigma$m at $\theta=\pi$. These
curves mimic the rg flow diagram shown in Fig. 1, and arrows on the
$\theta=0,\pi$ curves indicate rg flow from uv to ir. All curves agree in the
perturbative uv regime, while nonperturbative effects from the $\theta$ term
lead to divergent trajectories in the ir.

##  III Methods

In this work, we study the Hamiltonian defined in Eq. 3 by performing mc
sampling of the partition function $Z=\operatorname{Tr}e^{-\beta H}$ using a
worm algorithm [63, 64, 65] on a spacetime lattice. The dimensions
$L_{X}\times L_{Y}$ of the $2$-dimensional spatial lattices were varied in the
range $32\leq L_{X}\leq 1024$ and $L_{Y}=3,5,7$, with periodic boundary
conditions in $L_{X}$ and open boundary conditions in $L_{Y}$. The couplings
$J=|J^{\prime}|=1$ were held fixed and only the staggering $\gamma$ was varied
in the range $0\leq\gamma L_{Y}\leq 2$. The Heisenberg model with this level
of staggering is not frustrated and no sign problem occurs in the mc sampling.
In our computations, the imaginary-time extent $\beta$ was also discretized
into $L_{T}$ timesteps of size $\varepsilon$, such that $\beta=\varepsilon
L_{T}$. Strictly speaking, the Hamiltonian model of Eq. 3 is recovered only by
extrapolating to the $\varepsilon\to 0$ limit. Alternatively, one can develop
a cluster algorithm directly in the continuous time limit [66]. However, since
we are interested in studying the continuum limit of a relativistic field
theory, we perform mc computations at a fixed $\varepsilon=1.0$, which gives a
transfer matrix model with the same continuum limit. The physical extent
$\beta$ is chosen such that $\beta=L_{X}/c$, where $c$ is the speed of light
of this system.

For all combinations of the parameters we calculate the second-moment
correlation length $\xi_{2}(L_{X},g)$ from the spin-spin correlation function
$\langle\vec{S}_{x,y}\cdot\vec{S}_{x^{\prime},y^{\prime}}\rangle$. This long
distance length scale has been extensively studied at $\theta=0$ [52, 44] and
is easy to extract from lattice results. The calculation is then repeated with
doubled volume $2L_{X}$ but fixed bare couplings
$g_{\text{bare}}=(L_{Y},J^{\prime},J,\gamma)$. This macroscopic change in
scale $L_{X}\to 2L_{X}$ defines a discrete variant of the $\beta$-function,
known as the step-scaling function [42]

$\displaystyle F_{\xi}(z)$
$\displaystyle=\frac{\xi_{2}(2L_{X},g_{\text{bare}})}{\xi_{2}(L_{X},g_{\text{bare}})},$
$\displaystyle z=\xi_{2}(L_{X},g_{\text{bare}})/L_{X}.$ (5)

where $z$ defines a renormalized coupling. In the continuum limit
$\xi_{2}\to\infty$, at constant $z$, the step-scaling function $F_{\xi}(z)$
becomes a universal function, which uniquely characterizes the corresponding
qft.

##  IV Results

Figure 4 shows numerical results for the step-scaling function for the
$\text{O}(3)$ nl$\sigma$m at various $\theta$, computed using the qubit
Hamiltonian of Eq. 3. Results from both alternating (left panel) and columnar
(right panel) staggering configurations are shown. To guide the reader, we
show three continuous curves: the perturbative prediction (dotted line) and
nonperturbative mc results for $\theta=0$ (black dashed line), and
$\theta=\pi$ (black solid line). The perturbative curve is a two-loop
computation [44] valid in the uv ($z\gg 1$) and shows asymptotic freedom near
the uv fixed point $F(z)=2$ at $z\to\infty$. The $\theta=0$ curve was obtained
with the standard lattice action in Ref. [44], and shows the flow from the uv
to the trivial ir fixed point ($F_{\xi}(0)=1$ at $z=0$).

The $\theta=\pi$ curve (black solid line) is a polynomial fit in $z^{-2n}$ up
to order $n=5$ to our mc results with $\gamma L_{Y}=0$. This shows the rg flow
from the asymptotically-free uv fixed point at $z=\infty$ to the
$\text{SU}(2)_{1}$ wzw fixed point in the ir at $z=z^{*}$. We estimate the
location of the nontrivial ir fixed point to be $z^{\ast}\approx 0.28$ where
$F(z^{*})=2$, which is the discrete equivalent of a vanishing
$\beta$-function. We emphasize that the physics of _all_ scales, from
asymptotic freedom in the uv to the $\text{SU}(2)_{1}$ wzw theory in the ir is
reproduced by this model.

The remaining curves show new results for non-zero $\gamma
L_{Y}\sim|\theta-\pi|/\pi$. The $\theta=0,\pi$ curves form a lower and upper
bound on all step scaling curves $0\leq\theta\leq\pi$. All curves closely
follow the perturbative two-loop calculation (dotted line) at large $z$ down
to $z\approx 0.75$. At lower values of the renormalized coupling, non-
perturbative effects start to dominate, leading to divergent trajectories.

For small staggering $\gamma L_{Y}$ the curves closely track the $\theta=\pi$
curve. But since $\theta$ is a relevant perturbation about the wzw fixed
point, the rg trajectories cannot reach the nontrivial fixed point at
$z^{\ast}$ and ultimately have to flow away to the trivial fixed point at
$z=0$, consistent with the fact that these theories are massive. These
theories can be made to pass arbitrarily close to the $\text{SU}(2)_{1}$ wzw
fixed point by choosing smaller and smaller $\gamma L_{Y}$, without the need
for any fine tuning, exemplified by the curve $\gamma L_{Y}=0.01$ in Fig. 4.
This is the phenomenon of conformal walking, which is also exhibited by qcd-
like 4-dimensional non-Abelian gauge theories near the conformal window [67],
or technicolor extensions of the Standard model [47].

As the staggering $\gamma$ is increased further, the step scaling curves trace
out the entire area bounded by the two curves $\theta=0,\pi$, demonstrating
that all $\theta$ vacua are contained in this model. However, this only yields
a qualitative relationship between $\gamma L_{Y}$ and $\theta$. Semi-classical
results from large-$(SL_{Y})$ expansions [55, 68] suggest that the
relationship should be linear, $\theta=2\pi SL_{Y}(1+\gamma f(L_{Y}))$ with
$f(L_{Y})\to f(\infty)$ approaching a finite constant in the large-$L_{Y}$
limit. Numerically, we observe that a value of $\gamma L_{Y}=1.0$ ($\gamma
L_{Y}=0.25$) approximates the $\theta=0$ curve with the alternating (columnar)
staggering. Additional data also show a periodic reappearance of $\theta=\pi$
around values of $\gamma L_{Y}=2.0$ ($\gamma L_{Y}=0.5$). From this we
estimate the asymptotic values $f(\infty)\approx 1.0$ for alternating and
$f(\infty)\approx 0.25$ for columnar staggering. The small discontinuities of
the step-scaling curves between different values of $L_{Y}=3,5,7$ suggest that
the corrections to $f(L_{Y})$ at finite $L_{Y}$ are mild, especially
considering that similar values of the renormalized coupling
$z=\xi_{2}(L_{X},g)$ were obtained with drastically different lattice spacings
(usually $L_{X}=64$ for the smaller $L_{Y}$ compared to $L_{X}=1024$ for the
larger $L_{Y}$).

Similar results are also observed with even ladders, where the staggering
$\gamma$ is a perturbation about the $\theta=0$ theory. Preliminary results
from $L_{Y}=2,4,6$ with alternating and columnar staggerings suggest that
$\theta=\pi$ can also be obtained in this way. These results strongly motivate
a conjecture: the continuum limit ($L_{Y}\to\infty$, with $L_{Y}$ either odd
or even) for each fixed $\gamma L_{Y}$ is in fact a _unique_ qft corresponding
to the $1+1$-dimensional $\text{O}(3)$ nl$\sigma$m with a fixed $\theta$,

$\displaystyle\theta$ $\displaystyle\equiv\pi L_{Y}+\gamma f\quad\text{(mod
$2\pi$)},$ (6)

where $f$ is a non-universal constant which depends on the details of the
model such as choice of couplings, staggering configuration, and whether
$L_{Y}$ is odd or even. While we have provided strong evidence in favor of
this identification, there are many paths forward to establish this more
rigorously. For example, odd and even $L_{Y}$ could be used to self-validate
this conjecture, by showing that their step-scaling functions agree by
appropriately tuning $\gamma L_{Y}$. Further, a comparison with the approach
of Ref. [46] using topological lattice actions would be very illuminating. In
that approach, $\theta$ appears as a manifestly topological parameter and thus
does not require an empirical identification like Eq. 6. It would also be
interesting to connect with analytical results on the $\theta$ vacua based on
semiclassical instanton methods.

##  V Conclusions

In this work, we have shown how to implement the $1+1$-dimensional
$\text{O}(3)$ nl$\sigma$m at arbitrary $\theta$ using qubit degrees of
freedom. While the motivation behind this work is the quantum simulation of
$\theta$ vacua on near-term quantum hardware, interestingly, this result also
advances lattice computations of qfts using _classical_ mc methods. On the
classical side, it provides the first sign-problem-free mc algorithm for
arbitrary $\theta$. Our numerical results, obtained with an efficient worm
algorithm, indicate that the entire range of $\theta$ vacua is contained in
this model, and we conjecture a simple prescription of how the continuum limit
can be reached and the physics at all scales can be studied.

This construction enables real-time simulation of $\theta$ vacua in the
$\text{O}(3)$ nl$\sigma$m on near-term quantum hardware. These theories can be
regularized at any lattice spacing through an embedding into a $2$-dimensional
square lattice of qubits with nearest-neighbor Heisenberg-type interactions.
The alternating staggering is a prime candidate for an analog quantum
simulation platform like ultracold atoms, with uniform pairwise interactions,
and couplings that can be arranged through the trapping pattern shown in Fig.
3. On digital quantum hardware like superconducting qubits or trapped ions,
either staggering can be implemented using standard Suzuki-Trotter
decompositions. Interestingly, the limit $L_{X}\gg L_{Y}$ is also amenable to
dmrg-type algorithms on tensor networks, which would be powerful complementary
approach to lattice mc and quantum simulation going forward.

In lattice field theory, the $1+1$-dimensional $\text{O}(3)$ nl$\sigma$m has
been long considered an ideal testbed for static properties of qcd, exhibiting
many of its features, including asymptotic-freedom and $\theta$ vacua. Even
more possibilities open up once we have access to realtime dynamics using
quantum platforms. For instance, accessing nontrivial $\theta$ would allow the
study of inelastic scattering processes in an asymptotically-free theory,
which would have been impossible in the integrable $\theta=0,\pi$ theories.

Formulating qfts using qubits can yield unexpected advantages. For the
$\text{O}(3)$ nl$\sigma$m, this approach has the rather remarkable feature
that it completely circumvents a sign problem present in conventional lattice
formulations of the $\theta$ term, and is amenable to efficient cluster
algorithms. Extension to the entire family of $\text{CP}(N-1)$ models is
straightforward. This is encouraging on the path forward towards studying qcd
with novel classical and quantum algorithms. Our results demonstrate that
there is no fundamental obstruction to studying $\theta$ vacua with discrete
degrees of freedom, but whether such ideas might one day even help with the
sign problems in qcd remains to be seen.

## Acknowledgments

We learned about D-theory from Shailesh Chandrasekharan and Uwe-Jens Wiese and
are grateful to them for many enlightening conversations over the years. We
thank Martin Savage for inspiring discussions and important feedback on the
manuscript. We would also like to acknowledge stimulating conversations with
Tanmoy Bhattacharya, Anthony Ciavarella, Mendel Nguyen and Mithat Ünsal on
related matters.

The material presented here was funded in part by the DOE QuantISED program
through the theory consortium “Intersections of QIS and Theoretical Particle
Physics” at Fermilab with Fermilab Subcontract No. 666484, in part by
Institute for Nuclear Theory with US Department of Energy Grant DE-
FG02-00ER41132, and in part by U.S. Department of Energy, Office of Science,
Office of Nuclear Physics, Inqubator for Quantum Simulation (IQuS) under Award
Number DOE (NP) Award DE-SC0020970.

## References

  * Ratti [2018] C. Ratti, Lattice QCD and heavy ion collisions: A review of recent progress, Reports on Progress in Physics 81, 084301 (2018).
  * Davoudi  _et al._ [2021] Z. Davoudi, W. Detmold, P. Shanahan, K. Orginos, A. Parreño, M. J. Savage, and M. L. Wagman, Nuclear matrix elements from lattice QCD for electroweak and beyond-Standard-Model processes, Physics Reports Nuclear Matrix Elements from Lattice QCD for Electroweak and beyond–Standard-Model Processes, 900, 1 (2021).
  * Hansen and Sharpe [2019] M. T. Hansen and S. R. Sharpe, Lattice QCD and Three-Particle Decays of Resonances, Annual Review of Nuclear and Particle Science 69, 65 (2019).
  * Constantinou  _et al._ [2021] M. Constantinou, A. Courtoy, M. A. Ebert, M. Engelhardt, T. Giani, T. Hobbs, T.-J. Hou, A. Kusina, K. Kutak, J. Liang, H.-W. Lin, K.-F. Liu, S. Liuti, C. Mezrag, P. Nadolsky, E. R. Nocera, F. Olness, J.-W. Qiu, M. Radici, A. Radyushkin, A. Rajan, T. Rogers, J. Rojo, G. Schierholz, C. P. Yuan, J.-H. Zhang, and R. Zhang, Parton distributions and lattice-QCD calculations: Toward 3D structure, Progress in Particle and Nuclear Physics 121, 103908 (2021).
  * Hsu and Reeb [2010] S. D. Hsu and D. Reeb, On the sign problem in dense qcd, International Journal of Modern Physics A 25, 53 (2010).
  * Goy  _et al._ [2017] V. Goy, V. Bornyakov, D. Boyda, A. Molochkov, A. Nakamura, A. Nikolaev, and V. Zakharov, Sign problem in finite density lattice qcd, Progress of Theoretical and Experimental Physics 2017, 031D01 (2017).
  * Jordan  _et al._ [2011] S. P. Jordan, K. S. Lee, and J. Preskill, Quantum computation of scattering in scalar quantum field theories, arXiv preprint arXiv:1112.4833 (2011).
  * Jordan  _et al._ [2012] S. P. Jordan, K. S. Lee, and J. Preskill, Quantum algorithms for quantum field theories, Science 336, 1130 (2012).
  * Yeter-Aydeniz  _et al._ [2019] K. Yeter-Aydeniz, E. F. Dumitrescu, A. J. McCaskey, R. S. Bennink, R. C. Pooser, and G. Siopsis, Scalar quantum field theories as a benchmark for near-term quantum computers, Physical Review A 99, 032306 (2019).
  * Klco and Savage [2019] N. Klco and M. J. Savage, Digitization of scalar fields for quantum computing, Physical Review A 99, 052335 (2019).
  * Chandrasekharan  _et al._ [2002] S. Chandrasekharan, B. Scarlet, and U. J. Wiese, From spin ladders to the 2D O(3) model at non-zero density, Computer Physics Communications Proceedings of the Europhysics Conference on Computational Physics Computational Modeling and Simulation of Complex Systems, 147, 388 (2002).
  * Brower  _et al._ [2004] R. Brower, S. Chandrasekharan, S. Riederer, and U. J. Wiese, D-theory: Field quantization by dimensional reduction of discrete variables, Nuclear Physics B 693, 149 (2004).
  * Beard  _et al._ [2006] B. B. Beard, M. Pepe, S. Riederer, and U.-J. Wiese, Efficient Cluster Algorithm for CP(N-1) Models, Computer Physics Communications 175, 629 (2006).
  * Laflamme  _et al._ [2016] C. Laflamme, W. Evans, M. Dalmonte, U. Gerber, H. Mejía-Díaz, W. Bietenholz, U. J. Wiese, and P. Zoller, CP(N-1) quantum field theories with alkaline-earth atoms in optical lattices, Annals of Physics 370, 117 (2016).
  * Evans  _et al._ [2018] W. Evans, U. Gerber, M. Hornung, and U. J. Wiese, SU(3) quantum spin ladders as a regularization of the CP(2) model at non-zero density: From classical to quantum simulation, Annals of Physics 398, 94 (2018).
  * Bruckmann  _et al._ [2019] F. Bruckmann, K. Jansen, and S. Kühn, O(3) nonlinear sigma model in $1+1$ dimensions with matrix product states, Physical Review D 99, 074501 (2019).
  * Singh and Chandrasekharan [2019] H. Singh and S. Chandrasekharan, A qubit regularization of the ${O(3)}$ sigma model, Phys.Rev. D100, 054505 (2019).
  * Bhattacharya  _et al._ [2021] T. Bhattacharya, A. J. Buser, S. Chandrasekharan, R. Gupta, and H. Singh, Qubit regularization of asymptotic freedom, Physical Review Letters 126, 172001 (2021).
  * Singh [2019] H. Singh, Qubit $o(n)$ nonlinear sigma models, arXiv:1911.12353 [hep-lat, physics:quant-ph] (2019).
  * Singh [2022] H. Singh, Large-charge conformal dimensions at the $O(N)$ Wilson-Fisher fixed point, arXiv:2203.00059 [hep-lat] (2022).
  * Chandrasekharan and Wiese [1997] S. Chandrasekharan and U. J. Wiese, Quantum link models: A discrete approach to gauge theories, Nuclear Physics B 492, 455 (1997).
  * Brower  _et al._ [1999] R. Brower, S. Chandrasekharan, and U. J. Wiese, QCD as a quantum link model, Phys.Rev. D60, 094502 (1999).
  * Raychowdhury and Stryker [2020] I. Raychowdhury and J. R. Stryker, Solving gauss’s law on digital quantum computers with loop-string-hadron digitization, Physical Review Research 2, 033039 (2020).
  * Anishetty  _et al._ [2009] R. Anishetty, M. Mathur, and I. Raychowdhury, Prepotential formulation of su (3) lattice gauge theory, Journal of Physics A: Mathematical and Theoretical 43, 035403 (2009).
  * Banerjee  _et al._ [2013] D. Banerjee, M. Bögli, M. Dalmonte, E. Rico, P. Stebler, U.-J. Wiese, and P. Zoller, Atomic quantum simulation of u (n) and su (n) non-abelian lattice gauge theories, Physical review letters 110, 125303 (2013).
  * Zohar  _et al._ [2015] E. Zohar, J. I. Cirac, and B. Reznik, Quantum simulations of lattice gauge theories using ultracold atoms in optical lattices, Reports on Progress in Physics 79, 014401 (2015).
  * Banuls  _et al._ [2017] M. C. Banuls, K. Cichy, J. I. Cirac, K. Jansen, and S. Kühn, Efficient basis formulation for (1+ 1)-dimensional su (2) lattice gauge theory: spectral calculations with matrix product states, Physical Review X 7, 041046 (2017).
  * Muschik  _et al._ [2017] C. Muschik, M. Heyl, E. Martinez, T. Monz, P. Schindler, B. Vogell, M. Dalmonte, P. Hauke, R. Blatt, and P. Zoller, U (1) wilson lattice gauge theories in digital quantum simulators, New Journal of Physics 19, 103020 (2017).
  * Zache  _et al._ [2018] T. V. Zache, F. Hebenstreit, F. Jendrzejewski, M. Oberthaler, J. Berges, and P. Hauke, Quantum simulation of lattice gauge theories using wilson fermions, Quantum science and technology 3, 034010 (2018).
  * Alexandru  _et al._ [2019] A. Alexandru, P. F. Bedaque, S. Harmalkar, H. Lamm, S. Lawrence, N. C. Warrington, N. Collaboration, _et al._ , Gluon field digitization for quantum computers, Physical Review D 100, 114501 (2019).
  * Bender and Zohar [2020] J. Bender and E. Zohar, Gauge redundancy-free formulation of compact qed with dynamical matter for quantum and classical computations, Physical Review D 102, 114517 (2020).
  * Davoudi  _et al._ [2020] Z. Davoudi, M. Hafezi, C. Monroe, G. Pagano, A. Seif, and A. Shaw, Towards analog quantum simulations of lattice gauge theories with trapped ions, Physical Review Research 2, 023015 (2020).
  * Klco  _et al._ [2020] N. Klco, M. J. Savage, and J. R. Stryker, Su (2) non-abelian gauge field theory in one dimension on digital quantum computers, Physical Review D 101, 074512 (2020).
  * Shaw  _et al._ [2020] A. F. Shaw, P. Lougovski, J. R. Stryker, and N. Wiebe, Quantum algorithms for simulating the lattice schwinger model, Quantum 4, 306 (2020).
  * Kasper  _et al._ [2020] V. Kasper, T. V. Zache, F. Jendrzejewski, M. Lewenstein, and E. Zohar, Non-abelian gauge invariance from dynamical decoupling, arXiv preprint arXiv:2012.08620 (2020).
  * Buser  _et al._ [2021] A. J. Buser, H. Gharibyan, M. Hanada, M. Honda, and J. Liu, Quantum simulation of gauge theory via orbifold lattice, Journal of High Energy Physics 2021, 1 (2021).
  * Haase  _et al._ [2021] J. F. Haase, L. Dellantonio, A. Celi, D. Paulson, A. Kan, K. Jansen, and C. A. Muschik, A resource efficient approach for quantum and classical simulations of gauge theories in particle physics, Quantum 5, 393 (2021).
  * Zamolodchikov and Zamolodchikov [1979] A. B. Zamolodchikov and A. B. Zamolodchikov, Factorized S-matrices in two dimensions as the exact solutions of certain relativistic quantum field theory models, Annals of Physics 120, 253 (1979).
  * Balog and Hegedűs [2004] J. Balog and Á. Hegedűs, TBA equations for excited states in the O(3) andO(4) nonlinear -model, Journal of Physics A: Mathematical and General 37, 1881 (2004).
  * Balog and Hegedűs [2010] J. Balog and Á. Hegedűs, The finite size spectrum of the 2-dimensional O(3) nonlinear $\sigma$-model, Nuclear Physics B 829, 425 (2010).
  * Lüscher [1978] M. Lüscher, Quantum non-local charges and absence of particle production in the two-dimensional non-linear $\sigma$-model, Nuclear Physics B 135, 1 (1978).
  * Lüscher  _et al._ [1991] M. Lüscher, P. Weisz, and U. Wolff, A numerical method to compute the running coupling in asymptotically free theories, Nuclear Physics B 359, 221 (1991).
  * Bietenholz  _et al._ [1996] W. Bietenholz, A. Pochinsky, and U.-J. Wiese, Testing Haldane’s Conjecture in the O(3) Model by a Meron Cluster Simulation, Nuclear Physics B - Proceedings Supplements 47, 727 (1996).
  * Caracciolo  _et al._ [1995] S. Caracciolo, R. G. Edwards, A. Pelissetto, and A. D. Sokal, Asymptotic scaling in the two-dimensional O(3) $\sigma$ model at correlation length 1${0}^{5}$, Physical Review Letters 75, 1891 (1995).
  * Balog  _et al._ [2000] J. Balog, P. Forgacs, and L. Palla, A two-dimensional integrable axionic sigma-model and T-duality, arXiv:hep-th/0004180 10.1016/S0370-2693(00)00645-6 (2000).
  * Bogli  _et al._ [2012] M. Bogli, F. Niedermayer, M. Pepe, and U. J. Wiese, Non-trivial $\theta$-vacuum effects in the 2-d o(3) model, JHEP 04, 117.
  * de Forcrand  _et al._ [2012] P. de Forcrand, M. Pepe, and U. J. Wiese, Walking near a conformal fixed point: The 2-d $o(3)$ model at $\theta\mathbf{\approx}\pi$ as a test case, Physical Review D 86, 075006 (2012).
  * Wilson and Kogut [1974] K. G. Wilson and J. B. Kogut, The Renormalization group and the epsilon expansion, Phys.Rept. 12, 75 (1974).
  * Wilson [1983] K. G. Wilson, The renormalization group and critical phenomena, Reviews of Modern Physics 55, 583 (1983).
  * Kogut [1979] J. B. Kogut, An introduction to lattice gauge theory and spin systems, Reviews of Modern Physics 51, 659 (1979).
  * Beard  _et al._ [2005] B. B. Beard, M. Pepe, S. Riederer, and U.-J. Wiese, Study of $\mathrm{C}\mathrm{P}(n-1)$ $\theta$-vacua by cluster simulation of $\mathrm{S}\mathrm{U}(n)$ quantum spin ladders, Physical Review Letters 94, 010603 (2005).
  * Caracciolo  _et al._ [1993] S. Caracciolo, R. G. Edwards, A. Pelissetto, and A. D. Sokal, Wolff type embedding algorithms for general nonlinear sigma models, Nucl.Phys. B403, 475 (1993).
  * Shankar and Read [1990] R. Shankar and N. Read, The = $\pi$ nonlinear sigma model is massless, Nuclear Physics B 336, 457 (1990).
  * Affleck [1988] I. Affleck, Field Theory Methods and Quantum Critical Phenomena, in  _Les Houches Summer School in Theoretical Physics: Fields, Strings, Critical Phenomena_ (1988).
  * Martín-Delgado  _et al._ [1996] M. A. Martín-Delgado, R. Shankar, and G. Sierra, Phase Transitions in Staggered Spin Ladders, Physical Review Letters 77, 3443 (1996).
  * Sierra [1996] G. Sierra, The nonlinear sigma model and spin ladders, Journal of Physics A: Mathematical and General 29, 3299 (1996).
  * Haldane [1983a] F. D. M. Haldane, Continuum dynamics of the 1-D Heisenberg antiferromagnetic identification with the O(3) nonlinear sigma model, Phys.Lett. A93, 464 (1983a).
  * Haldane [1983b] F. D. M. Haldane, Nonlinear field theory of large spin Heisenberg antiferromagnets. Semiclassically quantized solitons of the one-dimensional easy Axis Neel state, Phys.Rev.Lett. 50, 1153 (1983b).
  * Affleck and Haldane [1987] I. Affleck and F. D. M. Haldane, Critical theory of quantum spin chains, Physical Review B 36, 5291 (1987).
  * Venuti  _et al._ [2005] L. C. Venuti, C. D. E. Boschi, E. Ercolessi, F. Ortolani, G. Morandi, S. Pasini, and M. Roncaglia, Particle content of the nonlinear sigma model with a $\theta$-term: A lattice model investigation, Journal of Statistical Mechanics: Theory and Experiment 2005, L02004 (2005).
  * Read and Sachdev [1989a] N. Read and S. Sachdev, Valence-bond and spin-Peierls ground states of low-dimensional quantum antiferromagnets, Physical Review Letters 62, 1694 (1989a).
  * Read and Sachdev [1989b] N. Read and S. Sachdev, Some features of the phase diagram of the square lattice SU(N) antiferromagnet, Nuclear Physics B 316, 609 (1989b).
  * Prokof’ev and Svistunov [2010] N. Prokof’ev and B. Svistunov, Worm Algorithm for Problems of Quantum and Classical Statistics, arXiv:0910.1393 [cond-mat, physics:hep-lat] (2010).
  * Prokof’ev and Svistunov [2001] N. Prokof’ev and B. Svistunov, Worm Algorithms for Classical Statistical Models, Phys.Rev.Lett. 87, 160601 (2001).
  * Wiese and Ying [1994] U. J. Wiese and H. P. Ying, A determination of the low energy parameters of the 2-d Heisenberg antiferromagnet, Zeitschrift für Physik B Condensed Matter 93, 147 (1994).
  * Beard and Wiese [1996] B. B. Beard and U.-J. Wiese, Simulations of Discrete Quantum Systems in Continuous Euclidean Time, Physical Review Letters 77, 5130 (1996).
  * Banks and Zaks [1982] T. Banks and A. Zaks, On the phase structure of vector-like gauge theories with massless fermions, Nuclear Physics B 196, 189 (1982).
  * Sierra [1997] G. Sierra, On the application of the nonlinear sigma model to spin chains and spin ladders, Lect.Notes Phys. 478, 137 (1997).