# Phase structure of the CP(1) model in the presence of a topological
$\theta$-term

Katsumasa Nakayama  NIC, DESY, Platanenallee 6, D-15738 Zeuthen, Germany  Lena
Funcke  Perimeter Institute for Theoretical Physics, 31 Caroline Street North,
Waterloo, ON N2L 2Y5, Canada  Karl Jansen  NIC, DESY, Platanenallee 6, D-15738
Zeuthen, Germany  Ying-Jer Kao  Department of Physics and Center for
Theoretical Physics, National Taiwan University, Taipei 10607, Taiwan  Physics
Division, National Center for Theoretical Science, National Taiwan University,
Taipei 10607, Taiwan  Stefan Kühn  Computation-Based Science and Technology
Research Center, The Cyprus Institute, 20 Kavafi Street, 2121 Nicosia, Cyprus

###### Abstract

We numerically study the phase structure of the CP(1) model in the presence of
a topological $\theta$-term, a regime afflicted by the sign problem for
conventional lattice Monte Carlo simulations. Using a bond-weighted Tensor
Renormalization Group method, we compute the free energy for inverse couplings
ranging from $0\leq\beta\leq 1.1$ and find a CP-violating, first-order phase
transition at $\theta=\pi$. In contrast to previous findings, our numerical
results provide no evidence for a critical coupling $\beta_{c}<1.1$ above
which a second-order phase transition emerges at $\theta=\pi$ and/or the
first-order transition line bifurcates at $\theta\neq\pi$. If such a critical
coupling exists, as suggested by Haldane’s conjecture, our study indicates
that is larger than $\beta_{c}>1.1$.

##  I Introduction

The CP($N-1$) models in 1+1 dimensions share many properties with QCD in 3+1
dimensions, among them confinement, asymptotic freedom, instantons, a $1/N$
expansion, a topological charge, and a $\theta$-term. Thus, they serve as
benchmark models for developing and testing both new numerical techniques and
new proposed solutions to open questions of QCD.

Both the CP($N-1$) models and QCD contain many non-perturbative phenomena,
which cannot be addressed with conventional lattice techniques. One particular
example is the above-mentioned $\theta$-term, which is the origin of the
strong CP problem in QCD Crewther  _et al._ (1979); Abel  _et al._ (2020). The
$\theta$-term gives rise to an imaginary contribution to the action in the
Euclidean formulation of lattice gauge theories. This complex action problem,
also called the sign problem Troyer and Wiese (2005), prevents the successful
application of Markov chain Monte Carlo (MCMC) methods for large values of the
topological vacuum angle $\theta$. To evade this problem, new numerical
techniques have to be developed.

One particularly promising approach are tensor network (TN) techniques, which
have been successfully employed to simulate lattice gauge theories in 1+1
dimensions Bañuls  _et al._ (2013a, b); Buyens  _et al._ (2014); Kühn  _et
al._ (2014); Buyens  _et al._ (2015, 2016a, 2016b); Bañuls  _et al._ (2017a,
b); Zapp and Orús (2017); Buyens  _et al._ (2017); Bañuls  _et al._ (2018);
Saito  _et al._ (2015a); Bañuls  _et al._ (2015); Saito  _et al._ (2015b);
Buyens  _et al._ (2016c); Bañuls  _et al._ (2016); Rico  _et al._ (2014);
Pichler  _et al._ (2016); Ercolessi  _et al._ (2018); Magnifico  _et al._
(2019a, b); Kühn  _et al._ (2015); Bañuls  _et al._ (2017c); Sala  _et al._
(2018, 2019); Silvi  _et al._ (2017, 2019); Shimizu and Kuramashi (2014a, b);
Kawauchi and Takeda (2016a, 2018); Byrnes  _et al._ (2002); Buyens  _et al._
(2017); Funcke  _et al._ (2020); Bruckmann  _et al._ (2019) and have recently
been applied to theories in 2+1 Tagliacozzo  _et al._ (2013, 2014); Zohar  _et
al._ (2015); Kuramashi and Yoshimura (2019); Felser  _et al._ (2020); Emonts
_et al._ (2020) and 3+1 Magnifico  _et al._ (2020) dimensions. In particular,
the $\theta$-dependence of lattice gauge theories has already been
successfully studied using various TN methods, i.e., the density matrix
renormalization group method Byrnes  _et al._ (2002), matrix product states
Buyens  _et al._ (2017); Funcke  _et al._ (2020), and the tensor
renormalization group (TRG) approach Kawauchi and Takeda (2016a, 2018);
Shimizu and Kuramashi (2014b). Specifically, TRG has been used to study the
$\theta$-dependence of the CP(1) model Kawauchi and Takeda (2016a, 2018) and
the Schwinger model Shimizu and Kuramashi (2014b).

The phase structure of different CP($N-1$) models has been intensively
investigated with different methods in the past Schierholz (1994); Plefka and
Samuel (1997a, b); Imachi  _et al._ (1999); Azcoiti  _et al._ (2007); Kawauchi
and Takeda (2016a, 2018). In particular, the phase diagram of the CP(1) model
with a $\theta$-term has been studied with a strong coupling analysis Plefka
and Samuel (1997b), Monte Carlo simulations Azcoiti  _et al._ (2007), and TRG
studies Kawauchi and Takeda (2016a, 2018). Here, it was shown that a first-
order phase transition occurs at small values of $\beta$ and $\theta=\pi$,
where the CP symmetry spontaneously breaks (see Fig. 1). For larger values of
$\beta$, corresponding to the weak-coupling regime, the numerical studies in
Refs. Azcoiti  _et al._ (2007); Kawauchi and Takeda (2016a, 2018) indicated
that the phase transition at $\theta=\pi$ becomes a second-order transition.
Since the CP(1) model is equivalent to the O(3) model Banerjee (1994), these
results are in agreement with Haldane’s conjecture Haldane (1983a, b) that the
O(3) model at weak coupling becomes gapless at $\theta=\pi$, corresponding to
a second-order phase transition. At the same time, it has also been suggested
that the first-order transition line bifurcates at $\theta\neq\pi$ for large
values $\beta$, thus exhibiting a non-trivial dependence on the coupling
Plefka and Samuel (1997b); Azcoiti  _et al._ (2007); Kawauchi and Takeda
(2016a, 2018). However, the exact value of the critical coupling, $\beta_{c}$,
remains unknown for both scenarios. While a strong coupling analysis suggests
a bifurcation at $\beta_{c}=0.56$ Plefka and Samuel (1997b), Monte Carlo
simulations yield a second-order transition beyond $\beta_{c}=0.5$ Azcoiti
_et al._ (2007), and the TRG analysis finds that the second-order phase
transition occurs around $\beta_{c}\approx 0.3-0.4$ Kawauchi and Takeda
(2016a, 2018).

In the present work, we reanalyze the phase diagram of the lattice CP(1) model
in the presence of a $\theta$-term using high-precision TRG simulations. We
provide a detailed analysis of the different errors contributing to the
current and previous numerical studies. Interestingly, we find no indication
of a bifurcation point or the onset of a second-order phase transition up to
$\beta=1.1$, which is much larger than suggested by the previous studies. If
the bifurcation point really exists, as suggested by Haldane’s conjecture, our
results indicate it occurs at $\beta_{c}>1.1$.

The paper is organized as follows. In Sec. II, we will describe the CP(1)
model and the numerical TRG methods. In Sec. III, we will present our
numerical results. In Sec. IV, we will summarize and discuss our results.

Figure 1: Sketch of the phase diagram of the lattice CP(1) model, where
$\beta=1/g^{2}$ is the inverse coupling and $\theta$ is the free parameter of
the topological $\theta$-term. The weak-coupling limit $\beta\to\infty$
corresponds to the continuum limit of the lattice model. For small values of
$\beta$, the model undergoes a first-order phase transition at $\theta=\pi$
(solid gray line), up to a certain critical coupling $\beta_{c}$ (yellow
shaded region). For $\beta>\beta_{c}$ (blue region), different predictions
exist and the phase structure is not entirely clear. There are hints towards a
second-order phase transition at $\theta=\pi$ (dotted brown line), in
agreement with Haldane’s conjecture, and towards a bifurcation of the critical
first-order line, such that the model undergoes two first-order transitions in
$\theta\neq\pi$ (dashed dark purple lines).

##  II Model and methods

In order to study the phase structure of the CP(1) model with the TRG Levin
and Nave (2007); Xie  _et al._ (2012); Kawauchi and Takeda (2016b), we use a
Euclidean time lattice formulation. The discretized action of the model is
given by

$S_{\theta}=-2\beta\sum_{x,\mu}\left[z^{*}_{x}z_{x+\hat{\mu}}U_{x,\mu}+z_{x}z^{*}_{x+\mu}U^{\dagger}_{x,\mu}\right]-i\frac{\theta}{2\pi}\sum_{x}q_{x},$
(1)

where $\beta=1/g^{2}$ is the inverse coupling constant and $\theta\in[0,2\pi]$
is the free parameter of the topological $\theta$-term. The two-component
complex scalar fields $z_{x}\in\mathds{C}^{2}$ reside on the sites
$x=(x_{1},x_{2})$ of a two-dimensional lattice and are fixed to unit length,
$z_{x}^{*}z_{x}=1$ $\forall x$. The gauge variables $U_{x,\mu}$ reside on the
links connecting the lattice sites $x$ and $x+\hat{\mu}$ along the direction
$\mu$ (see Fig. 2). They are related to the components $A_{x,\mu}$ of an
auxiliary vector field via $U_{x,\mu}=\exp(iA_{x,\mu})$. The topological
charge $q_{x}$ is defined as

$q_{x}=A_{x,1}-A_{x+\hat{1}-\hat{2},2}-A_{x-\hat{2},1}+A_{x-\hat{2},2}\mathrm{\
mod\ }2\pi.$ (2)

Figure 2: Illustration of the 1+1D square lattice showing the vertices (blue
spheres) and the directed links. The complex scalar fields $z_{x}$
(illustrated in black) are sitting on the vertices $x$ and the link variables
$U_{x,\mu}$ are sitting on the edges connecting sites $x$ and $x+\hat{mu}$ as
illustrated in green. The topological charge $q_{x}$ (orange) corresponds to a
product of the sum of the gauge fields around a plaquette.

In order to explore the phase diagram of the CP(1) model using TRG, we need to
obtain the tensor network representation of the partition function

$\displaystyle\mathcal{Z}_{\theta}=\int\prod_{x}\mathrm{d}z_{x}\prod_{y,\mu}\mathrm{d}A_{y,\mu}\
\exp(-S_{\theta}),$ (3)

where we have used a shorthand notation for the integration measure over the
complex scalar fields,

$\displaystyle\mathrm{d}z_{x}=\mathrm{d}^{2}z^{*}_{x}\ \mathrm{d}^{2}z_{x}\
\delta(|z_{x}|-1).$ (4)

To this end, we first derive an expansion for $\exp(-S_{\theta=0})$ which
factors into two parts: a part dependent on the complex scalar fields $z_{x}$
and a part depending on the vector fields $A_{x,\mu}$. Inserting this
expansion into Eq. (3) and taking the integral over the scalar fields, we
arrive at a tensor representation of the matter part. Subsequently, we derive
a corresponding representation for the $\theta$-term. The final tensors for
$\mathcal{Z}_{\theta}$ can then be obtained by multiplying both tensors for
the individual parts together and integrating over the gauge fields. To obtain
a tensor network representation for $\theta=0$, we use the character-like
expansion of the Boltzmann factor Plefka and Samuel (1997b)

$\displaystyle\begin{split}e^{-S_{\theta=0}}=\prod_{x,\mu}Z_{0}(\beta)\sum_{l,m=0}^{\infty}d_{l,m}\exp{[i(m-l)A_{x,{\mu}}]}\\\
\times h_{l,m}(\beta)f_{l,m}(z_{x},z_{x+\hat{\mu}}),\end{split}$ (5)

where $d_{l,m}$ are the expansion dimensions, $h_{l,m}$ are the expansion
coefficients, $f_{l,m}$ are the expansion characters, and $Z_{0}(\beta)$ is a
normalization factor. The non-negative integers $l$ and $m$ will become the
indices of the tensor representation after integrating out the degrees of
freedom $z_{x}$ and $A_{\mu}$. Note that due to the product in front of the
right-hand side of Eq. (5), there are four sets of indices $(l,m)$ associated
with each site. These correspond to the ingoing and outgoing index for every
direction (see Fig. 3(a)).

The expansion coefficients $h_{l,m}$ and the normalization factor
$Z_{0}(\beta)$ can be expressed through the modified Bessel functions of the
first kind, $I_{n}(x)$,

$\displaystyle h_{l,m}(\beta)$
$\displaystyle=\frac{I_{1+l+m}(4\beta)}{I_{1}(4\beta)},$ (6) $\displaystyle
Z_{0}(\beta)$ $\displaystyle=\frac{I_{1}(4\beta)}{2\beta}.$ (7)

As the explicit form of the dimensions $d_{l,m}$ and the characters $f_{l,m}$
is rather complicated, we refer to Refs. Plefka and Samuel (1997b); Kawauchi
and Takeda (2016b) for details and only quote their normalization conditions
here,

$\displaystyle\begin{split}\int\mathrm{d}z_{x}\
f_{l,m}(z_{x^{\prime}},z_{x})f^{*}_{l^{\prime},m^{\prime}}(z_{x^{\prime\prime}},z_{x})&=\frac{1}{d_{l,m}}f_{l,m}(z_{x^{\prime}},z_{x^{\prime\prime}}),\\\
f_{l,m}(z_{x},z_{x})&=d_{l,m},\end{split}$ (8)

In order to construct the tensor representation, we introduce an additional
decomposition of $f_{l,m}$ (see Fig. 3(b)),

$f_{l,m}(z_{x},z_{x^{\prime}})=\sum_{\\{a\\}}F_{l,m}^{\\{a\\}}(z_{x})\tilde{F}_{l,m}^{\\{a\\}}(z_{x^{\prime}}),$
(9)

where $\\{a\\}=a_{1},...,a_{l},a^{\prime}_{1},...,a^{\prime}_{m}$, and the
indices $a_{i}$ ($a_{i}^{\prime}$) range over $1,2$.

With this expansion, we can divide $f_{l,m}(z_{x},z_{x^{\prime}})$ into a
$z_{x}$-dependent part $F_{l,m}^{\\{a\\}}(z_{x})$ and a
$z_{x^{\prime}}$-dependent part $\tilde{F}_{l,m}^{\\{a\\}}(z_{x^{\prime}})$,
such that we can subsequently integrate out the complex scalar fields $z_{x}$
site by site. Note that the $\theta$-term does not depend on $z_{x}$, which
implies that only the first term in Eq. (1) is relevant for integrating out
$z_{x}$. Finally, we define the tensor (see Fig. 3(b)) by the integral

$\displaystyle\begin{aligned}
W^{x}&{}_{(l_{s},m_{s},\\{a\\})(l_{t},m_{t},\\{b\\})(l_{u},m_{u},\\{c\\})(l_{v},m_{v},\\{d\\})}\\\
\equiv&\,\sqrt{d_{l_{s},m_{s}}d_{l_{t},m_{t}}d_{l_{u},m_{u}}d_{l_{v},m_{v}}}\\\
&\times\,\sqrt{h_{l_{s},m_{s}}(\beta)h_{l_{t},m_{t}}(\beta)h_{l_{u},m_{u}}(\beta)h_{l_{v},m_{v}}(\beta)}\\\
&\times\,\int\mathrm{d}z_{x}\
\tilde{F}^{a}_{l_{s},m_{s}}(z_{x}){F}^{b}_{l_{t},m_{t}}(z_{x})\tilde{F}^{c}_{l_{u},m_{u}}(z_{x}){F}^{d}_{l_{v},m_{v}}(z_{x}).\end{aligned}$
(10)

The tensor $W_{x}$ represents a rank-4 tensor, whose bonds are given by the
four multi indices $(l_{s},m_{s},\\{a\\})$, $(l_{t},m_{t},\\{b\\})$,
$(l_{u},m_{u},\\{c\\})$, and $(l_{v},m_{v},\\{d\\})$, which allows for
representing $\exp(-S_{\theta=0})$ as a tensor network up to the factor
$e^{i(m-l)A_{x,\mu}}$ depending on the gauge field.

Next, we derive a tensor representation of the topological $\theta$-term
starting again from its character expansion, which reads

$e^{i\frac{\theta}{2\pi}q_{x}}=\sum_{n_{p}\in\mathbb{Z}}e^{in_{p}q_{x}}G_{n_{p}}(\theta).$
(11)

In the expression above, $n_{p}$ is an integer, and we have defined

$G_{n_{p}}(\theta)=\frac{2\,\mathrm{sin}\frac{\theta+2\pi
n_{p}}{2}}{\theta+2\pi n_{p}}.$ (12)

In order to take the integral over the gauge field $A_{x,\mu}$, we consider
the following integrals of $A_{x,1}$ and $A_{x-\hat{2},2}$,

$\displaystyle\begin{split}\int^{\pi}_{-\pi}\mathrm{d}A_{x,1}e^{i(m_{t}-l_{t})A_{x,1}}e^{i(t_{p}-v_{p})A_{x,1}}&=\delta^{l_{t}-m_{t}}_{t_{p}-v_{p}},\\\
\int^{\pi}_{-\pi}\mathrm{d}A_{x-\hat{2},2}e^{i(m_{u}-l_{u})A_{x-\hat{2},2}}e^{i(s_{p}-t_{p})A_{x-\hat{2},2}}&=\delta^{l_{u}-m_{u}}_{s_{p}-t_{p}}.\end{split}$
(13)

Similarly, the integral over $A_{x,2}$ yields
$\delta^{l_{t}-m_{t}}_{u_{p}-v_{p}}$. Using
$\delta_{ab}\delta_{bc}=\delta_{ac}$, we finally get the following tensor
representation,

$\displaystyle T^{x}_{stuv}$ $\displaystyle\equiv$ $\displaystyle
T^{x}_{(l_{s},m_{s},\\{a\\},s_{p})(l_{t},m_{t},\\{b\\},t_{p})(l_{u},m_{u},\\{c\\},u_{p})(l_{v},m_{v},\\{d\\},v_{p})}$
(14) $\displaystyle\equiv$ $\displaystyle
W^{x}_{(l_{s},m_{s},\\{a\\})(l_{t},m_{t},\\{b\\})(l_{u},m_{u},\\{c\\})(l_{v},m_{v},\\{d\\})}$
$\displaystyle\times$
$\displaystyle\delta^{l_{t}-m_{t}}_{t_{p}-v_{p}}\delta^{l_{u}-m_{u}}_{s_{p}-t_{p}}G_{t_{p}}(\theta)\delta_{t_{p}u_{p}}.$

The partition function of the system can now be represented as the contraction
of the rank-4 tensors $T_{stuv}^{x}$ sitting at the different lattice sites.
Compared to $W^{x}$, the multi-index for each bond of the tensor $T^{x}$
comprises an additional index resulting from the character expansion of the
topological term (see also Fig. 3(c)).

Since the indices $l$, $m$, and $n_{p}$ in the character expansions in Eq. (1)
and Eq. (11) range over a countable infinite set, the bonds of the tensor
$T^{x}$ are infinite dimensional. Hence, for numerical calculations we need to
truncate the bonds to a finite dimension. To this end, we limit the number of
terms that we keep in the character expansions. For the indices $l$, $m$ we
define a cutoff $k_{\mathrm{max}}$ and restrict ourselves to values fulfilling
$l+m\leq k_{\mathrm{max}}$. This results in a bond size of
$\chi_{\beta}=1+k_{\mathrm{max}}2^{k_{\mathrm{max}}+1}$ for the tensor
representing the partition function without the $\theta$-term. Similarly, we
only keep a subset of $\chi_{\theta}$ terms in the expansion for the
topological term in Eq. (11). Since our numerical analysis primarily focuses
on the region $\theta=\pi$, we choose the terms with the largest absolute
value of $G_{n_{p}}(\pi)$. For example, for $\chi_{\theta}=2$, we take into
account $n_{p}=0,-1$ because the absolute value of $|G_{0,-1}(\pi)|={2}/{\pi}$
is larger than all other values of $|G_{n_{p}}(\pi)|$. The size of the
truncated tensor $T^{x}$ is thus given by $\chi_{\beta}\times\chi_{\theta}$.

These cutoffs in the tensor $T^{x}$ introduce systematic errors. Since
$h_{l,m}$ is given by a ratio of the modified Bessel functions (see Eq. (6)),
the cutoff at $k_{\mathrm{max}}$ truncates only an exponentially small
contribution to the tensor $T_{stuv}^{x}$. The truncation in the expansion for
the $\theta$-term leads to an error of the order
$\mathcal{O}(1/\chi_{\theta})$, as shown in Eq. (12). We will estimate these
systematic errors by using different $k_{\mathrm{max}}$ and $\chi_{\theta}$
our numerical calculations.

Figure 3: (a) Illustration of the function $f_{l,m}(z_{x},z_{x+\hat{\mu}})$
(indicated by the shaded region behind the gauge links) connecting the complex
scalar field at site $x$ (blue sphere) to its neighboring lattice sites (gray
spheres). Note that each site has two ingoing and two outgoing links, hence
there are four pairs of integer indices $(l_{i},m_{i})$, $i=s,t,u,v$ at each
site. (b) Illustration of the situation after using Eq. (8) to break the
function $f_{l,m}(z_{x},z_{x+\hat{\mu}})$ into a contraction of two tensors
between neighboring sites. The colored part belongs to site $x$, the gray
parts belong to the neighboring sites. The gray square illustrates the pieces
which form the tensor $W^{x}$ in Eq. (10) upon integrating over the scalar
field. Panel (c) shows the final rank-4 tensor $T^{x}_{stuv}$ sitting at
vertex $x$. The bonds of the tensor are represented as legs sticking out of
the vertex. Each bond corresponds to a multi-index corresponding to a pair of
indices $(l_{i},m_{i})$ resulting from the character expansion (solid black
and red lines), two sets of indices $\\{a\\}$ from the decomposition of the
function $f_{l_{i},m_{i}}$ in Eq. (9) (dashed black and red lines), and an
index from the expansion of topological $\theta$-term (black dotted line).

In addition, the TRG algorithm introduces another systematic error. Starting
from the initial, finite dimensional tensor $T^{x}$, the coarse-gaining step
of the algorithm would produce a tensor which has the same structure as the
original one, but with a bond size corresponding to the square of the one of
the original tensor. Hence, for the computation to be sustainable, one has to
truncate the tensor during the coarse-graining step and limit the bond size to
a maximum value $D$. We will estimate this additional systematic error by
performing calculations for various values of the bond dimensions $D$ for a
fixed set of parameters. In particular, $D$ should be chosen larger than the
size of the initial tensor $\chi_{\beta}\times\chi_{\theta}$ to capture the
relevant physics of the model.

##  III Results

In order to explore the phase structure of the model, we compute the partition
function using TRG methods, which in turn allows us to obtain the free energy
density

$\displaystyle F(\beta,\theta)=-\frac{1}{\beta V}\log\mathcal{Z}_{\theta},$
(15)

where $V$ is the dimensionless volume of the lattice. Moreover, we can also
compute the specific heat from the partition function using the relation

$\displaystyle
C(\beta,\theta)=\frac{\beta^{2}}{V}\frac{\partial^{2}}{\partial\beta^{2}}\log\mathcal{Z}_{\theta}.$
(16)

In particular, studying the free energy will allow us to determine the phase
structure of the model and to detect a possible phase transition as
discontinuities in the derivatives of the free energy.

As a first step, we focus on the case $\theta=0$ to explore systematic errors
due to the truncation in $k_{\mathrm{max}}$ over a large range $\beta$. This
allows us to identify the range of couplings in which our TRG computations
yield reliable results. As a second step, we turn to the model in the presence
of a topological $\theta$-term. Using again the TRG approach, we study the
model in a region around $\theta=\pi$, where we carefully estimate our
systematic errors due to the finite values of
$(D,k_{\mathrm{max}},\chi_{\theta},V)$.

###  III.1 CP(1) model without $\theta$-term

To begin with, we benchmark the TRG method for $\theta=0$ and explore the
effect of systematic errors due to the truncation of the expansion in Eq. (5).
To this end, we calculate the free energy and specific heat of the CP(1) model
using the Anisotropic TRG (ATRG) method with randomized singular value
decomposition (SVD) Morita  _et al._ (2018); Adachi  _et al._ (2020). In order
to avoid errors due to the randomized SVD, we set the oversampling parameter
to the value of $D$, meaning that we compute $2D$ singular values at each
coarse-graining step and select the $D$ dominant ones from those. This ATRG
method is faster than the usual TRG algorithm, while the precision is
comparable to the one of TRG for same $D$. Since we choose the oversampling
parameter large enough, the systematic error from the randomized SVD is
negligible.

Figure 4 shows our results for the free energy and the specific heat for
$D=80$, $k_{\mathrm{max}}=2,4$, and a volume of $V=2^{40}$.

Figure 4: ATRG results for the free energy $F(\beta,\theta=0)$ (main plot) and
the specific heat $C(\beta,\theta=0)$ (inset) as a function of the inverse
coupling $\beta=1/g^{2}$, with the volume $V=2^{40}$, bond size $D=80$, and
cutoff $k_{\mathrm{max}}=2$ (blue dots) and $4$ (orange triangles).

Focusing on the free energy, we see that the systematic error due to the
finite value of $k_{\mathrm{max}}$ is negligible in the range of small values
of the inverse coupling. In particular, for $\beta\leq 1.5$ there is
essentially no difference between results with $k_{\mathrm{max}}=2$ and $4$.
Only as we enter the region of larger $\beta$, the curves start to deviate
slightly while the qualitative behavior of the free energy is still the same
for both values of $k_{\mathrm{max}}$.

Looking at the specific heat, in Fig. 4, we see that truncation effects due to
finite $k_{\mathrm{max}}$ are more pronounced. This is not too surprising, as
the specific heat is obtained via Eq. (16) by approximating the second
derivative with finite differences,

$\displaystyle\begin{split}C(&\beta,\theta=0)\approx\frac{\beta^{2}}{V}\frac{1}{12(\Delta\beta)^{2}}\times\\\
&\Bigl{[}16\bigl{(}\log\mathcal{Z}_{\theta}|_{\beta+\Delta\beta}\bigr{.}-2\log\mathcal{Z}_{\theta}|_{\beta}\bigl{.}+\log\mathcal{Z}_{\theta}|_{\beta-\Delta\beta}\bigr{)}\\\
&-\bigl{(}\log\mathcal{Z}_{\theta}|_{\beta+2\Delta\beta}-2\log\mathcal{Z}_{\theta}|_{\beta}\bigr{.}\Bigl{.}\bigl{.}+\log\mathcal{Z}_{\theta}|_{\beta-2\Delta\beta}\bigr{)}\Bigr{]},\end{split}$
(17)

where we choose $\Delta\beta=0.5$ for the data shown in the inset of Fig. 4.
Thus, the systematic errors in $\log\mathcal{Z}_{\theta}$ are enhanced by a
factor of $(\Delta\beta)^{2}$ leading to significantly larger deviations in
the regime of larger $\beta$. Focusing on smaller values of the inverse
coupling, our data show that truncation effects are nevertheless small despite
this effect, and we can reliably determine the specific heat with
$k_{\mathrm{max}}=2$ in regimes $\beta\leq 1.1$. Hence, for all the following
we will solely focus on that region and use $k_{\mathrm{max}}\geq 2$, for
which systematic errors due to the truncation of the expansion in Eq. (5) are
very small.

We note that previous work using TRG Kawauchi and Takeda (2018) observed large
fluctuations in the specific heat in the regime of large $\beta$. In contrast,
our ATRG results are stable even for large values of $\beta$. While Ref.
Kawauchi and Takeda (2018) used a bond dimension $D=21$, our calculations rely
on a larger bond dimension $D=80$. Thus, we suspect that these fluctuations
were caused by truncation effects due to the small bond dimension.

###  III.2 CP(1) model with $\theta$-term

In order to study the CP(1) model in the presence of a $\theta$-term, we use
the bond-weighted TRG method Adachi  _et al._ (2020) in conjunction with
regular SVD to truncate the tensor at each coarse-graining step. Compared to
the original TRG method, bond-weighted TRG introduces bond weights on the
edges of the tensor network, allowing for more precise results within the same
calculation time.

Figure 5 shows our results for the free energy as a function of $\theta$ for
the volume $V=2^{24}$, bond dimension $D=80$, $k_{\mathrm{max}}=2$, a
truncation of the character expansion of the topological term
$\chi_{\theta}=2$, and three values of the inverse coupling
$\beta=0.1,0.6,1.1$. Focusing first on $\beta=0.1$ in Fig. 5(a), we observe a
clear cusp at $\theta=\pi$. Hence, the first derivative of
$F(\beta=0.1,\theta)$ with respect to $\theta$ is discontinuous at
$\theta=\pi$, indicating that the model undergoes a first-order phase
transition at this point. This observation is consistent with previous
findings of a first-order transition for $\beta=0.1$ Plefka and Samuel
(1997b); Azcoiti  _et al._ (2007); Kawauchi and Takeda (2016a, 2018).

For larger inverse couplings of $\beta=0.6$ and $\beta=1.1$, the absolute
values of the free energy change, as Fig. 5(b) and Fig. 5(c) reveal,
respectively. However, for both $\beta=0.6$ and $\beta=1.1$, we still find the
same cusp structure as for the smaller inverse coupling of $\beta=0.1$. Thus,
even for these large values of $\beta$, the transition is still of first
order. These results disagree with previous findings Plefka and Samuel
(1997b); Azcoiti  _et al._ (2007); Kawauchi and Takeda (2016a, 2018) that for
$\beta_{c}\gtrsim 0.4-0.56$, the transition at $\theta=\pi$ becomes of second
order and the first-order transition line bifurcates at $\theta\neq\pi$.

Figure 5: TRG results for the free energy as a function of $\theta$, for the
volume $V=2^{24}$, the bond dimension $D=80$, and the inverse couplings (a)
$\beta=0.1$, (b) $\beta=0.6$, and (c) $\beta=1.1$. The vertical dashed gray
line indicates $\theta=\pi$, at which we expect the phase transition to occur.
Figure 6: TRG results for the systematic errors of computing the free energy
as a function of $\theta$, for the inverse couplings of (a) $\beta=0.1$, (b)
$\beta=0.6$, and (c) $\beta=1.1$. The blue dots in panels (a) and (b) show
$\Delta F(\beta,\theta)$ for $D=112$, $k_{\mathrm{max}}=2$ and
$\chi_{\theta}=2$, the orange triangles for $D=144$, $k_{\mathrm{max}}=3$ and
$\chi_{\theta}=2$, and the green squares for $D=144$, $k_{\mathrm{max}}=2$ and
$\chi_{\theta}=4$. For panel (c), the markers encode the same values of
$k_{\mathrm{max}}$ and $\chi_{\theta}$, but we keep a fixed value of $D=112$
in all cases.

In order to ensure that the discrepancies for the phase structure obtained
from our numerical results and previous studies are not caused by systematic
errors due to the choice of $k_{\mathrm{max}}$, $D$, and $\chi_{\theta}$, we
carefully examine the effect of these parameters. Figure 6 shows the
difference in the results for the free energy,

$\displaystyle\Delta
F(\beta,\theta)=\left|F(\beta,\theta)|_{(D,k_{\mathrm{max}},\chi_{\theta})}-F(\beta,\theta)|_{(80,2,2)}\right|,$
(18)

between $(D,k_{\mathrm{max}},\chi_{\theta})=(80,2,2)$ and various larger
choices of $k_{\mathrm{max}}$, $D$, and $\chi_{\theta}$ for a fixed volume of
$V=2^{24}$. Note that changing $\chi_{\theta}$ and $k_{\mathrm{max}}$ affects
the size of the initial tensor as discussed above. While for
$(k_{\mathrm{max}},\chi_{\theta})=(2,2)$ the bond size of the initial tensor
is $34$, for $(2,4)$ and $(3,2)$ it grows to $68$ and $98$ respectively. Since
$D$ should be chosen well above the initial tensor size, we do not only
increase $k_{\mathrm{max}}$ or $\chi_{\theta}$, but simultaneously also the
value of $D$.

Figure 6(a) shows our results for the smallest value of the inverse coupling
corresponding to $\beta=0.1$. In particular, we observe that increasing the
bond dimension while keeping $(k_{\mathrm{max}},\chi_{\theta})=(2,2)$ has an
extremely small effect on the results, with an absolute change in the range
below $10^{-14}$. Similarly, increasing $k_{\mathrm{max}}$ does not lead to
noticeable deviations and the absolute change is in the same range. Increasing
the value of $\chi_{\theta}$ to 4 yields a slightly larger deviation, which is
still very small and in the range of $10^{-10}$. All in all, when we compare
these absolute differences to the values of the free energy, which are on the
order of $10^{-1}$ (see Fig. 5 (a)), the systematic errors due to finite $D$,
$k_{\mathrm{max}}$ and $\chi_{\theta}$ are negligible throughout the entire
range of $\theta$ we study.

For the larger inverse couplings of $\beta=0.6$ and $\beta=1.1$ in Fig. 6(b)
and Fig. 6(c), respectively, we observe a qualitatively similar picture. While
the absolute deviation is larger than for the case of $\beta=0.1$, these
effects are well below the percentile range even for our largest value of
$\beta$. Just as before, increasing $D$ and $k_{\mathrm{max}}$ has the least
impact, whereas increasing $\chi_{\theta}$ to $4$ now clearly has a slightly
larger effect on the results. In relation to the values of the free energy
shown in Fig. 5(b) and Fig. 5(c), which are on the order of $10^{-1}$ and 3,
respectively, the systematic effects due to finite
$(D,k_{\mathrm{max}},\chi_{\theta})$ are again negligible.

Figure 7: TRG results for the free energy as a function of $\theta$ in a
region very close to $\theta=\pi$ (indicated by the vertical gray dashed
line), for the bond dimension $D=80$ and the inverse couplings (a)
$\beta=0.1$, (b) $\beta=0.6$, and (c) $\beta=1.1$. The different markers
correspond to the volumes $V=2^{12}$ (blue dots), $2^{14}$ (orange triangles),
and $2^{24}$ (green squares.)

To further study the possible origin of the discrepancies between our and
previous numerical results, we now investigate the cusp in the free energy as
a function of the volume $V$. To this end, we plot a close-up around the
region of $\theta=\pi$ for various volumes and various values of the inverse
coupling in Fig. 7. Looking at our results for $\beta=0.1$ (see Fig. 7(a)), we
observe that for volumes below $2^{14}$, the free energy is smooth at
$\theta=\pi$. Thus, the first derivative does not show any discontinuity for
small volumes. Only when increasing $V$ to larger values, the cusp at
$\theta=\pi$ eventually emerges, signaling the onset of the first-order
transition directly in the free energy. For larger values of the inverse
coupling, corresponding to $\beta=0.6$ and $1.1$, we see a qualitatively
similar behavior, as Figs. 7(b) and 7(c) reveal. Interestingly, the finite-
volume effects seem to be slightly less pronounced for larger values of
$\beta$, as a comparison between Figs. 7(a) and 7(c) shows. For $\beta=1.1$,
we observe that the data for $V=2^{14}$ already has a clear signature of a
cusp, and increasing the volume only leads to slightly sharper peak.

In summary, the analysis of the systematic errors from finite
$(V,D,k_{\mathrm{max}},\chi_{\theta})$ further supports the reliability of our
findings throughout the entire range of $\beta$ that we study. We conclude
that the phase transition at $\theta=\pi$ is still of first order up to an
inverse coupling of $\beta=1.1$.

##  IV Discussion and Conclusions

In this paper, we numerically studied the phase diagram of the lattice CP(1)
model with a $\theta$-term. We observe a CP-violating, first-order phase
transition at $\theta=\pi$ and small $\beta$, in agreement with previous
findings Plefka and Samuel (1997b); Azcoiti  _et al._ (2007); Kawauchi and
Takeda (2016a, 2018). However, our results do not confirm previous indications
for a bifurcation in the first-order transition line at some critical value of
$\beta$, ranging from $\beta_{c}=0.56$ Plefka and Samuel (1997b) to
$\beta_{c}=0.5$ Azcoiti  _et al._ (2007) to $\beta_{c}=0.4$ Kawauchi and
Takeda (2016a, 2018). Instead, we find that the first-order transition at
$\theta=\pi$ persists up to $\beta=1.1$ without any bifurcation, and we
observe no indication that the phase transition at $\theta=\pi$ becomes of
second order for $\beta>\beta_{c}$. The existence of a second-order transition
at $\theta=\pi$ and $\beta>\beta_{c}$ was suggested by the previous studies
based on a strong coupling analysis Plefka and Samuel (1997b), Monte Carlo
simulations Azcoiti  _et al._ (2007), and first TRG studies Kawauchi and
Takeda (2016a, 2018).

Compared to the strong coupling analysis of the CP(1) model in Ref. Plefka and
Samuel (1997b), our method does not require an expansion in the inverse
coupling $\beta=1/g^{2}$, but instead performs truncations and character
expansions to construct the tensor representation. Thus, while for small
$\beta$ the results in Plefka and Samuel (1997b) may be more accurate than
ours, the strong coupling analysis becomes more challenging for larger values
of $\beta$. Since our TRG results clearly show that the effects due to the
truncations of the character expansions are very small throughout the entire
range of $\beta$ we study, our results suggest that the strong coupling
expansion does not correctly capture the physics of the model up to
$\beta=1.1$.

Similarly, the Monte Carlo studies in Ref. Azcoiti  _et al._ (2007) suffer
from the sign problem, which makes it difficult to reliably examine the phase
diagram of the CP(1) model for large values of $\theta$. In particular, the
Monte Carlo approach becomes more challenging when approaching the phase
transition at $\theta=\pi$, whereas our TRG calculation only shows small
systematic errors, which are essentially independent of the value of $\theta$
for the parameter range we study.

Compared to previous TRG-based studies of the model, our approach to determine
the phase structure is more direct. In particular, Refs. Kawauchi and Takeda
(2016a, 2018) did not study the free energy of the model at large volumes, but
rather explored the free energy at small volumes and used finite-size scaling
behavior of the topological susceptibility as an indicator to determine the
order of the phase transition. In contrast, our computations are performed at
a much larger volume, for which the results are well converged, which allows
us to study the behavior of the free energy without having to rely on finite-
size scaling. Indeed, as we have shown in Fig. 7, a clear signature of the
first-order transition only occurs at volumes of $V>2^{14}$. We note that our
bond dimension of $D=80$ is only slightly larger than the bond dimension
$D=68$ used in Ref. Kawauchi and Takeda (2018), but increasing the bond
dimension only has a minor effect on our results, see Fig. 6.

To summarize, our results demonstrate that the phase transition at
$\theta=\pi$ is of first order up to $\beta\leq 1.1$, and no second-order
transition occurs. The discovery of such a second-order transition would be
crucial for confirming Haldane’s conjecture Haldane (1983a, b) that the O(3)
model becomes gapless at $\theta=\pi$ and weak coupling. As our study reveals,
the possibility of a second-order transition at $\theta=\pi$ and a bifurcation
of the first-order transition line at $\theta\neq\pi$ is only given for
$\beta_{c}>1.1$. Compared to previous studies, our method does not suffer from
the sign problem as Monte Carlo simulations, and our results have small
systematic errors of $\lesssim 10^{-3}$ due to larger volumes than in earlier
TRG studies. Thus, our study reveals that testing Haldane’s conjecture in the
weak-coupling regime requires both major numerical efforts and detailed
systematic error analysis. In the future, we will extend our study to larger
values of $\beta>1.1$. While this task is more challenging due to enhanced
truncation effects in the weak-coupling regime, our numerical make us positive
that this parameter range can be reliably accessed with TRG methods.

###### Acknowledgements.

Research at Perimeter Institute is supported in part by the Government of
Canada through the Department of Innovation, Science and Industry Canada and
by the Province of Ontario through the Ministry of Colleges and Universities.
S.K. acknowledges financial support from the Cyprus Research and Innovation
Foundation under project “Future-proofing Scientific Applications for the
Supercomputers of Tomorrow(FAST)”, contract no. COMPLEMENTARY/0916/0048.

## References

  * Crewther  _et al._ (1979) R. J. Crewther, P. Di Vecchia, G. Veneziano, and E. Witten, Chiral Estimate of the Electric Dipole Moment of the Neutron in Quantum Chromodynamics, Phys. Lett. B 88, 123 (1979), [Erratum: Phys.Lett.B 91, 487 (1980)].
  * Abel  _et al._ (2020) C. Abel _et al._ (nEDM), Measurement of the permanent electric dipole moment of the neutron, Phys. Rev. Lett. 124, 081803 (2020).
  * Troyer and Wiese (2005) M. Troyer and U.-J. Wiese, Computational Complexity and Fundamental Limitations to Fermionic Quantum Monte Carlo Simulations, Phys. Rev. Lett. 94, 170201 (2005).
  * Bañuls  _et al._ (2013a) M. C. Bañuls, K. Cichy, K. Jansen, and J. I. Cirac, The mass spectrum of the Schwinger model with Matrix Product States, J. High Energy Phys. 2013 (11), 158.
  * Bañuls  _et al._ (2013b) M. C. Bañuls, K. Cichy, J. I. Cirac, K. Jansen, and H. Saito, Matrix Product States for Lattice Field Theories, Proc. Sci. LATTICE2013 332, (2013b).
  * Buyens  _et al._ (2014) B. Buyens, J. Haegeman, K. Van Acoleyen, H. Verschelde, and F. Verstraete, Matrix Product States for Gauge Field Theories, Phys. Rev. Lett. 113, 091601 (2014).
  * Kühn  _et al._ (2014) S. Kühn, J. I. Cirac, and M. C. Bañuls, Quantum simulation of the Schwinger model: A study of feasibility, Phys. Rev. A 90, 042305 (2014).
  * Buyens  _et al._ (2015) B. Buyens, J. Haegeman, K. Van Acoleyen, and F. Verstraete, Matrix product states for Hamiltonian lattice gauge theories, Proc. Sci. LATTICE2014 308, (2015).
  * Buyens  _et al._ (2016a) B. Buyens, J. Haegeman, F. Verstraete, and K. Van Acoleyen, Tensor networks for gauge field theories, Proc. Sci., LATTICE2015 280, (2016a).
  * Buyens  _et al._ (2016b) B. Buyens, J. Haegeman, H. Verschelde, F. Verstraete, and K. Van Acoleyen, Confinement and String Breaking for ${Q}{E}{D}_{2}$ in the Hamiltonian Picture, Phys. Rev. X 6, 041040 (2016b).
  * Bañuls  _et al._ (2017a) M. C. Bañuls, K. Cichy, J. I. Cirac, K. Jansen, and S. Kühn, Density Induced Phase Transitions in the Schwinger Model: A Study with Matrix Product States, Phys. Rev. Lett. 118, 071601 (2017a).
  * Bañuls  _et al._ (2017b) M. C. Bañuls, K. Cichy, J. I. Cirac, K. Jansen, S. Kühn, and H. Saito, The multi-flavor Schwinger model with chemical potential - Overcoming the sign problem with Matrix Product States, Proc. Sci. LATTICE2016 316, (2017b).
  * Zapp and Orús (2017) K. Zapp and R. Orús, Tensor network simulation of QED on infinite lattices: Learning from (1+1) d, and prospects for (2+1) d, Phys. Rev. D 95, 114508 (2017).
  * Buyens  _et al._ (2017) B. Buyens, S. Montangero, J. Haegeman, F. Verstraete, and K. Van Acoleyen, Finite-representation approximation of lattice gauge theories at the continuum limit with tensor networks, Phys. Rev. D 95, 094509 (2017).
  * Bañuls  _et al._ (2018) M. C. Bañuls, K. Cichy, J. I. Cirac, K. Jansen, and S. Kühn, Tensor Networks and their use for Lattice Gauge Theories, Proc. Sci., LATTICE2018 334, (2018).
  * Saito  _et al._ (2015a) H. Saito, M. C. Bañuls, K. Cichy, J. I. Cirac, and K. Jansen, The temperature dependence of the chiral condensate in the Schwinger model with Matrix Product States, Proc. Sci. LATTICE2014 302, (2015a).
  * Bañuls  _et al._ (2015) M. C. Bañuls, K. Cichy, J. I. Cirac, K. Jansen, and H. Saito, Thermal evolution of the Schwinger model with Matrix Product Operators, Phys. Rev. D 92, 034519 (2015).
  * Saito  _et al._ (2015b) H. Saito, M. C. Bañuls, K. Cichy, J. I. Cirac, and K. Jansen, Thermal evolution of the one-flavour Schwinger model using Matrix Product States, Proc. Sci. LATTICE2015 283, (2015b).
  * Buyens  _et al._ (2016c) B. Buyens, F. Verstraete, and K. Van Acoleyen, Hamiltonian simulation of the Schwinger model at finite temperature, Phys. Rev. D 94, 085018 (2016c).
  * Bañuls  _et al._ (2016) M. C. Bañuls, K. Cichy, K. Jansen, and H. Saito, Chiral condensate in the Schwinger model with matrix product operators, Phys. Rev. D 93, 094512 (2016).
  * Rico  _et al._ (2014) E. Rico, T. Pichler, M. Dalmonte, P. Zoller, and S. Montangero, Tensor Networks for Lattice Gauge Theories and Atomic Quantum Simulation, Phys. Rev. Lett. 112, 201601 (2014).
  * Pichler  _et al._ (2016) T. Pichler, M. Dalmonte, E. Rico, P. Zoller, and S. Montangero, Real-Time Dynamics in U(1) Lattice Gauge Theories with Tensor Networks, Phys. Rev. X 6, 011023 (2016).
  * Ercolessi  _et al._ (2018) E. Ercolessi, P. Facchi, G. Magnifico, S. Pascazio, and F. V. Pepe, Phase transitions in ${Z}_{n}$ gauge models: Towards quantum simulations of the Schwinger-Weyl QED, Phys. Rev. D 98, 074503 (2018).
  * Magnifico  _et al._ (2019a) G. Magnifico, D. Vodola, E. Ercolessi, S. P. Kumar, M. Müller, and A. Bermudez, Symmetry-protected topological phases in lattice gauge theories: Topological ${\mathrm{QED}}_{2}$, Phys. Rev. D 99, 014503 (2019a).
  * Magnifico  _et al._ (2019b) G. Magnifico, D. Vodola, E. Ercolessi, S. P. Kumar, M. Müller, and A. Bermudez, $\mathbb{Z}_{N}$ gauge theories coupled to topological fermions: QED2 with a quantum-mechanical $\theta$ angle, Phys. Rev. B 100, 115152 (2019b).
  * Kühn  _et al._ (2015) S. Kühn, E. Zohar, J. Cirac, and M. C. Bañuls, Non-Abelian string breaking phenomena with matrix product states, J. High Energy Phys. 2015 (7), 130.
  * Bañuls  _et al._ (2017c) M. C. Bañuls, 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, Phys. Rev. X 7, 041046 (2017c).
  * Sala  _et al._ (2018) P. Sala, T. Shi, S. Kühn, M. C. Bañuls, E. Demler, and J. I. Cirac, Variational study of U(1) and SU(2) lattice gauge theories with Gaussian states in 1+1 dimensions, Phys. Rev. D 98, 034505 (2018).
  * Sala  _et al._ (2019) P. Sala, T. Shi, S. Kühn, M. C. Bañuls, E. Demler, and J. I. Cirac, Gaussian states for the variational study of (1+1)-dimensional lattice gauge models, Proc. Sci., LATTICE2018 334, (2019).
  * Silvi  _et al._ (2017) P. Silvi, E. Rico, M. Dalmonte, F. Tschirsich, and S. Montangero, Finite-density phase diagram of a (1+1)-d non-abelian lattice gauge theory with tensor networks, Quantum 1, 9 (2017).
  * Silvi  _et al._ (2019) P. Silvi, Y. Sauer, F. Tschirsich, and S. Montangero, Tensor network simulation of an su(3) lattice gauge theory in 1d, Phys. Rev. D 100, 074512 (2019).
  * Shimizu and Kuramashi (2014a) Y. Shimizu and Y. Kuramashi, Grassmann tensor renormalization group approach to one-flavor lattice Schwinger model, Phys. Rev. D 90, 014508 (2014a).
  * Shimizu and Kuramashi (2014b) Y. Shimizu and Y. Kuramashi, Critical behavior of the lattice Schwinger model with a topological term at $\theta=\pi$ using the Grassmann tensor renormalization group, Phys. Rev. D 90, 074503 (2014b).
  * Kawauchi and Takeda (2016a) H. Kawauchi and S. Takeda, Phase structure analysis of CP(N-1) model using Tensor renormalization group, PoS LATTICE2016, 322 (2016a).
  * Kawauchi and Takeda (2018) H. Kawauchi and S. Takeda, Loop-TNR analysis of CP(1) model with theta term, EPJ Web Conf. 175, 11015 (2018).
  * Byrnes  _et al._ (2002) T. M. R. Byrnes, P. Sriganesh, R. J. Bursill, and C. J. Hamer, Density matrix renormalization group approach to the massive Schwinger model, Phys. Rev. D 66, 013002 (2002).
  * Funcke  _et al._ (2020) L. Funcke, K. Jansen, and S. Kühn, Topological vacuum structure of the Schwinger model with matrix product states, Phys. Rev. D 101, 054507 (2020).
  * 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, Phys. Rev. D 99, 074501 (2019).
  * Tagliacozzo  _et al._ (2013) L. Tagliacozzo, A. Celi, A. Zamora, and M. Lewenstein, Optical Abelian Lattice Gauge Theories, Annals Phys. 330, 160 (2013).
  * Tagliacozzo  _et al._ (2014) L. Tagliacozzo, A. Celi, and M. Lewenstein, Tensor Networks for Lattice Gauge Theories with continuous groups, Phys. Rev. X 4, 041024 (2014).
  * Zohar  _et al._ (2015) E. Zohar, M. Burrello, T. Wahl, and J. I. Cirac, Fermionic Projected Entangled Pair States and Local U(1) Gauge Theories, Annals Phys. 363, 385 (2015).
  * Kuramashi and Yoshimura (2019) Y. Kuramashi and Y. Yoshimura, Three-dimensional finite temperature ${Z}_{2}$ gauge theory with tensor network scheme, J. High Energy Phys. 2019, 23.
  * Felser  _et al._ (2020) T. Felser, P. Silvi, M. Collura, and S. Montangero, Two-dimensional quantum-link lattice Quantum Electrodynamics at finite density, Phys. Rev. X 10, 041040 (2020).
  * Emonts  _et al._ (2020) P. Emonts, M. C. Bañuls, J. I. Cirac, and E. Zohar, Variational Monte Carlo simulation with tensor networks of a pure $\mathbb{Z}_{3}$ gauge theory in (2+1)d, Phys. Rev. D 102, 074501 (2020).
  * Magnifico  _et al._ (2020) G. Magnifico, T. Felser, P. Silvi, and S. Montangero, Lattice Quantum Electrodynamics in (3+1)-dimensions at finite density with Tensor Networks, (2020), arXiv:2011.10658 [hep-lat] .
  * Schierholz (1994) G. Schierholz, Towards a dynamical solution of the strong CP problem, Nucl. Phys. B Proc. Suppl. 37, 203 (1994).
  * Plefka and Samuel (1997a) J. C. Plefka and S. Samuel, Monte Carlo studies of two-dimensional systems with a theta term, Phys. Rev. D 56, 44 (1997a).
  * Plefka and Samuel (1997b) J. C. Plefka and S. Samuel, A Strong coupling analysis of the lattice CPN-1 models in the presence of a $\theta$ term, Phys. Rev. D 55, 3966 (1997b).
  * Imachi  _et al._ (1999) M. Imachi, S. Kanou, and H. Yoneyama, Two-dimensional CP2 model with $\theta$-term and topological charge distributions, Prog. Theor. Phys. 102, 653 (1999).
  * Azcoiti  _et al._ (2007) V. Azcoiti, G. Di Carlo, and A. Galante, Critical Behaviour of CP1 at $\theta=\pi$, Haldane’s Conjecture, and the Relevant Universality Class, Phys. Rev. Lett. 98, 257203 (2007).
  * Banerjee (1994) R. Banerjee, Quantum equivalence of O(3) nonlinear $\sigma$ model and the CP1 model: A Gauge independent Hamiltonian approach, Phys. Rev. D 49, 2133 (1994).
  * Haldane (1983a) F. D. M. Haldane, Continuum dynamics of the 1-D Heisenberg antiferromagnetic identification with the O(3) nonlinear sigma model, Phys. Lett. A 93, 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).
  * Levin and Nave (2007) M. Levin and C. P. Nave, Tensor renormalization group approach to 2D classical lattice models, Phys. Rev. Lett. 99, 120601 (2007).
  * Xie  _et al._ (2012) Z. Y. Xie, J. Chen, M. P. Qin, J. W. Zhu, L. P. Yang, and T. Xiang, Coarse-graining renormalization by higher-order singular value decomposition, Phys. Rev. B (2012).
  * Kawauchi and Takeda (2016b) H. Kawauchi and S. Takeda, Tensor renormalization group analysis of CP(N-1) model, Phys. Rev. D 93, 114503 (2016b).
  * Morita  _et al._ (2018) S. Morita, R. Igarashi, H.-H. Zhao, and N. Kawashima, Tensor renormalization group with randomized singular value decomposition, Phys. Rev. E 97, 033310 (2018).
  * Adachi  _et al._ (2020) D. Adachi, T. Okubo, and S. Todo, Anisotropic tensor renormalization group, Phys. Rev. B 102, 054432 (2020).
  * Adachi  _et al._ (2020) D. Adachi, T. Okubo, and S. Todo, Bond-weighted Tensor Renormalization Group, arXiv:2011.01679 (2020).