#  Simulation of many body localization and time crystals in two dimensions
with the neighborhood tensor update

Jacek Dziarmaga<EMAIL_ADDRESS>Jagiellonian University, Institute
of Theoretical Physics, ul. Łojasiewicza 11, 30-348 Kraków, Poland

###### Abstract

The Heisenberg antiferromagnet with discrete disorder on an infinite square
lattice is evolved in time from an initial Néel state. The simulation is
performed with the neighborhood tensor update (NTU) algorithm for an infinite
projected entangled pair state (iPEPS) [Phys. Rev. B 104, 094411 (2021)].
Ancillary spins are used to average over $2$ or $5$ discrete values of
disorder. With a bond dimension up to $20$, evolution times are long enough to
identify a many body localized regime for a strong enough disorder.
Furthermore, the same Hamiltonian is subject to periodic spin flips.
Simulations of the Floquet dynamics show that it can sustain a time
crystalline stage for a strong enough disorder.

##  I Introduction

Weakly entangled quantum states can be represented efficiently by tensor
networks Verstraete  _et al._ (2008); Orús (2014), including the one-
dimensional (1D) matrix product state (MPS) Fannes  _et al._ (1992), its two-
dimensional (2D) generalization known as a projected entangled pair state
(PEPS) Nishio  _et al._ (2004); Verstraete and Cirac (2004a), or a multi-scale
entanglement renormalization ansatz Vidal (2007, 2008); Evenbly and Vidal
(2014a, b). The MPS ansatz provides a compact representation of ground states
of 1D gapped local Hamiltonians Verstraete  _et al._ (2008); Hastings (2007);
Schuch  _et al._ (2008) and purifications of their thermal states Barthel
(2017). It is also the ansatz underlying the density matrix renormalization
group (DMRG) White (1992, 1993); Schollwöck (2005); Schöllwock (2011).
Analogously, the 2D PEPS is expected to represent ground states of 2D gapped
local Hamiltonians Verstraete  _et al._ (2008); Orús (2014) and their thermal
states Wolf  _et al._ (2008); Molnar  _et al._ (2015), though representability
of area-law states, in general, was shown to have its limitations Ge and
Eisert (2016). As a variational ansatz tensor networks do not suffer from the
sign problem plaguing the quantum Monte Carlo. Consequently, they can deal
with fermionic systems Corboz  _et al._ (2010a); Pineda  _et al._ (2010);
Corboz and Vidal (2009); Barthel  _et al._ (2009); Gu  _et al._ (2010), as was
shown for both finite Kraus  _et al._ (2010) and infinite PEPS (iPEPS) Corboz
_et al._ (2010b, 2011).

The PEPS was originally proposed as an ansatz for ground states of finite
systems Verstraete and Cirac (2004b); Murg  _et al._ (2007), generalizing
earlier attempts to construct trial wave-functions for specific models Nishio
_et al._ (2004). The subsequent development of efficient numerical methods for
an infinite PEPS (iPEPS) Jordan  _et al._ (2008); Jiang  _et al._ (2008); Gu
_et al._ (2008); Orús and Vidal (2009), shown in Fig. 1(a), promoted it as one
of the methods of choice for strongly correlated systems in 2D. Its power was
demonstrated, e.g., by a solution of the long-standing magnetization plateaus
problem in the highly frustrated compound
$\textrm{SrCu}_{2}(\textrm{BO}_{3})_{2}$ Matsuda  _et al._ (2013); Corboz and
Mila (2014), establishing the striped nature of the ground state of the doped
2D Hubbard model Zheng  _et al._ (2017), and new evidence supporting gapless
spin liquid in the kagome Heisenberg antiferromagnet Liao  _et al._ (2017).
Recent developments in iPEPS optimization Phien  _et al._ (2015); Corboz
(2016a); Vanderstraeten  _et al._ (2016), contraction Fishman  _et al._
(2018); Xie  _et al._ (2017), energy extrapolations Corboz (2016b), and
universality-class estimation Corboz  _et al._ (2018); Rader and Läuchli
(2018); Rams  _et al._ (2018) pave the way towards even more complicated
problems, including simulation of thermal states Czarnik  _et al._ (2012);
Czarnik and Dziarmaga (2014, 2015a); Czarnik  _et al._ (2016a); Czarnik and
Dziarmaga (2015b); Czarnik  _et al._ (2016b, 2017); Dai  _et al._ (2017);
Czarnik  _et al._ (2019a); Czarnik and Corboz (2019); Kshetrimayum  _et al._
(2019); Czarnik  _et al._ (2019b); Jiménez  _et al._ (2021, 2020); Poilblanc
_et al._ (2020); Czarnik  _et al._ (2021); Poilblanc  _et al._ (2021), mixed
states of open systems Kshetrimayum  _et al._ (2017); Czarnik  _et al._
(2019a); Mc Keever and Szymańska (2021), excited states Vanderstraeten  _et
al._ (2015); Ponsioen and Corboz (2020), or real-time evolution Czarnik  _et
al._ (2019a); Hubig and Cirac (2019); Hubig  _et al._ (2020); Abendschein and
Capponi (2008); Kshetrimayum  _et al._ (2020, 2021); Dziarmaga (2021); Schmitt
_et al._ (2021); Kaneko and Danshita (2021). In parallel with iPEPS, there is
continuous progress in simulating systems on cylinders of finite width using
DMRG. This numerically highly stable method that is now routinely used to
investigate 2D ground states Zheng  _et al._ (2017); Cincio and Vidal (2013)
was applied also to thermal states on a cylinder Bruognolo  _et al._ (2017);
Chen  _et al._ (2018a, 2019); Li  _et al._ (2019); Chen  _et al._ (2020).
However, the exponential growth of the bond dimension limits the cylinder’s
width to a few lattice sites. Among alternative approaches are direct
contraction and renormalization of a 3D tensor network representing a 2D
thermal density matrix Li  _et al._ (2011); Xie  _et al._ (2012); Ran  _et
al._ (2012, 2013, 2018); Peng  _et al._ (2017); Chen  _et al._ (2018b); Ran
_et al._ (2019) and a cluster expansion Vanhecke  _et al._ (2021).

In this paper we apply the recently defined neighborhood tensor update (NTU)
algorithm Dziarmaga (2021); Schmitt  _et al._ (2021) to simulate unitary time
evolution of a 2D many body localizing (MBL) system: the antiferromagnetic
spin-$1/2$ Heisenberg model with discrete disorder initialized in the Néel
state. The same simulatons were performed previously Hubig and Cirac (2019);
Kshetrimayum  _et al._ (2020) with the full update (FU) and simple update (SU)
algorithms. NTU was intended as a reasonable trade off between FU and SU which
is more accurate than SU but, at the same time, more efficient and stable than
FU. The efficiency allows to reach bond dimensions up to $20$ that in turn
permit longer evolution times. The longer evolution enables a more reliable
identification of a MBL regime for strong enough disorder. Encouraged by this
quantitative progress, we supplement the model with periodic spin flips. This
Floquet dynamics was previously simulated with SU Kshetrimayum  _et al._
(2021). Here we obtain long enough evolution times, converged in the bond
dimension, to identify time crystals for strong enough disorder.

This paper is organized as follows. In Sec. II we overview the NTU algorithm
whose more technical details can be found in App. A. The applications follow
in Sec. III and IV to, respectively, many body localization and time crystals.
We conclude in Sec. V.

Figure 1:  Essential NTU.  In (a) infinite PEPS with tensors $A$ (lighter
green) and $B$ (darker green) on the two sublattices. The red lines are
physical spin indices and the black lines are bond indices, with bond
dimension D, contracting NN sites. In one of Suzuki-Trotter steps a Trotter
gate is applied to every horizontal NN pair of $A$-$B$ tensors (but not to
horizontal $B$-$A$ pairs). The gate can be represented as a contraction of two
tensors by an index with dimension $r$. When the two tensors are absorbed into
tensors $A$ and $B$ the bond dimension between them increases from $D$ to
$r\times D$. In (b) the $A$-$B$ pair – with a Trotter gate applied to it – is
approximated by a pair of new tensors, $A^{\prime}$ (lighter blue) and
$B^{\prime}$ (darker blue), connected by an index with the original dimension
$D$. The new tensors are optimized to minimize the Frobenius norm of the
difference between the two networks in (b). The networks surround the
considered NN bond with its six NN tensors. They provide a minimal tensor
environment – necessary to make efficient use of limited $D$ – which can be
contracted exactly in an efficient way. After $A^{\prime}$ and $B^{\prime}$
are converged, they replace all tensors $A$ and $B$ in a new iPEPS. Then the
next Trotter gate can be applied. More details of the NTU algorithm can be
found in appendix A.

##  II Neighborhood tensor update

In the following NTU simulations we use the second order Suzuki-Trotter
decomposition of small time steps. An application of a two-site NN Trotter
gate is explained in a diagrammatic form in Fig. 1. NTU is intermediate
between the two most popular simulation schemes: the simple update (SU) and
full update (FU) Czarnik  _et al._ (2019a). In both after a Trotter gate is
applied to a pair of nearest neighbor (NN) sites a bond dimension of the index
between the sites is increased by a factor equal to the SVD rank of the gate.
In order to prevent its exponential growth with time the dimension is
truncated to a predefined value, $D$, in a way that minimizes error incurred
by the truncation. The two schemes differ by a measure of the error: FU takes
into account full infinite tensor environment while SU only the bonds adjacent
to the NN sites. The former is expected to perform better in case of long
range correlations while the latter is, at least formally, more efficient
thanks to its locality though it makes less efficient use of the limited bond
dimension than FU. In NTU the error measure is induced by the sites that are
NN to the Trotter gate. It was demonstrated Dziarmaga (2021) that in practice
it compromises the FU accuracy only a little for a price of small numerical
overhead over SU, hence it may turn out to be an optimal trade off for many
applications, especially when quantum correlations are not exceedingly long
like in e.g. Kibble-Zurek quenches in 2D Schmitt  _et al._ (2021). In this
sense it seems to be tailored for time evolution of MBL systems where the
localization is expected to limit the growth of correlations.

NTU can be placed in a broader context by noting that it is a special case of
a cluster update Wang and Verstraete (2011) where the size of the tensor
environment is a variable parameter interpolating between a local environment
in SU and an infinite one in FU. In case of ground state calculations the
cluster update was thoroughly investigated in Refs. Lubasch  _et al._ , 2014a,
b where an interplay between maximal achievable correlation length and the
cluster size was demonstrated. In NTU the cluster includes the neighboring
sites only, see Fig. 1, to allow the error measure to be calculated exactly
but with little numerical overhead over SU. The calculation involves only
tensor contractions that are fully parallelizable. Its exactness warrants the
error measure to be a manifestly Hermitian and non-negative quadratic form.
This property is essential for stability of NTU and makes it distinct from FU
where an approximate corner transfer matrix renormalization Orús and Vidal
(2009); Orús (2014) often breaks the Hermiticity and non-negativeness. In case
of longer correlations the small environment can, admittedly, make NTU
converge with the bond dimension more slowly than FU but this can be often
compensated by its better numerical efficiency and stability that allow NTU to
reach higher bond dimensions Dziarmaga (2021).

Figure 2:  Many body localization.  A spin imbalance in function of time for
different disorder strengths $h$ and ancillary dimensions $d_{A}$. Here the
many body localizing Hamiltonian (2) was initialized with the Néel state of
spins (3). For $d_{A}=2$ and $5$ we show $h=0,1,..,4$ and $h=0,1,...,8$,
respectively. The time step $dt=0.01$ and the evolution is terminated when the
Trotter gate error (11) exceeds $\delta=0.03$ for the first time.

##  III Many body localization

When excited the MBL systems evade thermalization but retaining memory of
their initial conditions Basko  _et al._ (2006); Abanin  _et al._ (2019). The
memory is sustained by localization that gives rise to quasi-local constants
of motion Pal and Huse (2010); Vosk  _et al._ (2015); Goihl  _et al._ (2018)
that allow the evolving state to equilibrate but to a non-thermal stationary
state. In 2D which is harder to investigate by any means — analytic, numerical
or experimental Bordia  _et al._ (2017) — the issue is far from settled. The
2D MBL was suggested to be unstable towards a crossover to ergodicity on very
long time scales De Roeck and Huveneers (2017). Several numerical approaches
were proposed to address the 2D problem Geißler and Pupillo (2020); Wahl  _et
al._ (2019); Kennes (2018); Hubig and Cirac (2019); De Tomasi  _et al._
(2019); Andraschko  _et al._ (2014); Paredes  _et al._ (2005); Kshetrimayum
_et al._ (2020). This paper follows the approach of simulating the time
evolution of a 2D MBL system with the iPEPS tensor network Kennes (2018);
Hubig and Cirac (2019); Andraschko  _et al._ (2014); Paredes  _et al._ (2005);
Kshetrimayum  _et al._ (2020) but it employs the NTU algorithm Dziarmaga
(2021). For a localizing system this algorithm is expected to offer an optimal
combination of stability and efficiency.

In order to be more specific, we consider the antiferromagnetic Heisenberg
model with disorder on an infinite square lattice:

$H=\sum_{\langle
j,j^{\prime}\rangle}\vec{S}_{j}\vec{S}_{j^{\prime}}+\sum_{j}h_{j}S^{z}_{j}.$
(1)

Here $\vec{S}_{j}=\frac{1}{2}\vec{\sigma}_{j}$ is a spin-$1/2$ operator at
site $j$. Every random $h_{j}$ is drawn independently from a set of $d_{A}$
evenly spaced discrete values in an interval $[-h,h]$. For instance,
$h_{j}\in\\{-h,h\\}$ for $d_{A}=2$, $h_{j}\in\\{-h,0,h\\}$ for $d_{A}=3$, etc.
The Hamiltonian is initialized with a Néel state, $|\text{N\'{e}el}\rangle$,
with spins pointing up/down along the $z$-axis on sublattice A/B of an
infinite checkerboard lattice. This initial product state is evolved in time.
Expectation values of spin operators at time $t$ are averaged over all
possible realisations of the disorder.

Here as in Refs. Hubig and Cirac, 2019; Andraschko  _et al._ , 2014; Paredes
_et al._ , 2005; Kshetrimayum  _et al._ , 2020 an equivalent formulation of
the problem is considered:

$H_{\rm MBL}=\sum_{\langle
j,j^{\prime}\rangle}\vec{S}_{j}\vec{S}_{j^{\prime}}+\frac{h}{S}\sum_{j}A^{z}_{j}S^{z}_{j},$
(2)

where $A^{z}_{j}$ is an ancilla spin-$S$ operator at site $j$ that takes
$d_{A}=2S+1$ evenly spaced discrete values in an interval $[-S,S]$. In this
formulation the initial state is replaced by an equivalent uniform product:

$\mbox{$|\psi(0)\rangle$}=\mbox{$|\text{N\'{e}el}\rangle$}\prod_{j}\mbox{$|+\rangle$}_{j},$
(3)

where $j$ runs over ancillas. Every ancillary spin-$S$ is prepared in a
superposition:

$\mbox{$|+\rangle$}_{j}=(2S+1)^{-1/2}\sum_{m_{j}=-S}^{S}\mbox{$|m_{j}\rangle$},$
(4)

with the same probability for all evenly spaced values of disorder
$h_{j}=m_{j}h/S\in[-h,h]$. The two formulations are equivalent in the sense
that a disorder-averaged expectation value of any operator $\cal{O}$ at time
$t$ in the former is equal to
$\mbox{$\langle\psi(t)|$}\cal{O}\mbox{$|\psi(t)\rangle$}$ in the latter.

The Hamiltonian (2) is the same as in Ref. Kshetrimayum  _et al._ , 2020 and
the same ancillary dimensions $d_{A}=2$ and $5$ are considered for a range of
disorder strengths. Unlike Ref. Kshetrimayum  _et al._ , 2020 NTU is employed
instead of SU and bond dimensions up to $D=20$ are reached instead of $D=4,5$.
As a result NTU evolution times are long enough to reach some semi-
quantitative conclusions without any extrapolation in time.

Figure 3:  Time crystallization.  A stroboscopic time evolution of
magnetization on sublattice $A$ in function of number of cycles for two
combinations $(d_{A},h)$ of ancillary dimension and disorder strength,
$(2,50)$ and $(5,100)$, and two values of the angle deficit parameter:
$\epsilon=0$ and $0.5$. Here the Floquet Hamiltonian (5) with period $T=0.1$
was initialized with a Néel state of spins along $S^{Z}$ (3). The time step
was $dt=0.001$ and the evolution was terminated when the Trotter gate error
$\delta$, defined in (11), exceeded for the first time the value in the
bottom-right corner of each panel. In the top and bottom-right panels this
value is chosen such that the results for $D=18,20$ appear converged until the
termination. In the bottom-left case there is still some dependence on $D$
with a tendency to improve the time crystal with increasing $D$.

Figure 2 shows time evolution of spin imbalance between sublattices $A$ and
$B$,
$\mbox{$\langle\psi(t)|$}S^{z}_{A}\mbox{$|\psi(t)\rangle$}-\mbox{$\langle\psi(t)|$}S^{z}_{B}\mbox{$|\psi(t)\rangle$}$,
in function of evolution time $t$ for several values of disorder strength $h$.
In particular both panels show the evolution without disorder, $h=0$, when the
imbalance crosses from positive to negative. Not quite surprisingly, in the
absence of MBL this simulation becomes poorely converged in $D$ when the
system begins to thermalize and growing entanglement becomes difficult to
accommodate by the tensor network. The same comment, though to a lesser
extent, applies to a finite but weak disorder with $h=1$. For higher disorder
strengths the difference between $d_{A}=2$ and $d_{A}=5$ appears to be rather
quantitative. In the former case, as far as the limited evolution time permits
to conclude, there seems to be a crossover to MBL between $h=1$ and $h=2$. For
$h\geq 2$ the spin imbalance at first goes down from the initial $1$ before it
bounces up. The memory of the initial spin imbalance between the sublattices
is not lost. In the case of $d_{A}=5$ a similar bounce up is observed, within
available evolution time, only for the highest disorder strengths $h=7,8$. For
$h=2,...,6$ the NTU evolution terminates before bouncing up. It remains an
open question whether the spin imbalance is bouncing up at later times or the
evolution is heading straight towards thermalization and how precisely it
depends on $h$. For the available evolution times we can crudely estimate the
crossover between the (hypothetical) thermalization and MBL to take place
between $h=1$ and $h=7$.

It is worth emphasizing that, thanks to NTU (instead of SU) and bond
dimensions up to $D=20$ (instead of $D=4,5$), our ($D$-converged) evolution
times are typically twice as long as in Ref. Kshetrimayum  _et al._ , 2020.
This allows us to identify the MBL regime not only for $d_{A}=5$ but also for
the two-level disorder.

In the next section the same Hamiltonian (2) is promoted to a Floquet
Hamiltonian by the action of time-periodic spin flips.

Figure 4:  Time crystalization failure.  Two generic examples of stroboscopic
time evolution where the disorder is too weak to sustain a time crystal (TC)
that extends beyond the achievable simulation time. In the left panel TC
cannot be identified at all while in the right one there seems to be a
transient TC stage that begins to disappear after around $40$ Floquet cycles.
Here the Floquet Hamiltonian (5) with period $T=0.1$ was initialized with the
Néel state of spins (3). The time step was $dt=0.001$ and the evolution was
terminated when the Trotter gate error $\delta$, defined in (11), exceeded for
the first time the value in the bottom-right corner of each panel.

##  IV Time crystals

Time crystals were at first envisioned by Wilczek as spontaneous time
translation symmetry breaking in a ground state of a quantum Hamiltonian
Wilczek (2012) but this original idea was proved wrong by a no go theorem
Watanabe and Oshikawa (2015); Kozin and Kyriienko (2019). More recently
spontaneous discrete time translation symmetry breaking was demonstrated Sacha
(2015); Khemani  _et al._ (2016); Else  _et al._ (2016); Yao  _et al._ (2017)
in Floquet systems that are periodic in time. These ideas inspired experiments
that were performed in a variety of physical platforms Zhang  _et al._ (2017);
Choi  _et al._ (2017); Rovny  _et al._ (2018a); Pal  _et al._ (2018); Rovny
_et al._ (2018b); Smits  _et al._ (2018); Liao  _et al._ (2019); Smits  _et
al._ (2020). The fast growing field is already subject to several reviews
Sacha and Zakrzewski (2017); Khemani  _et al._ (2019); Else  _et al._ (2020).
A few mechanisms were proposed how a Floquet system can avoid heating. They
include the many body localization Ponte  _et al._ (2015); Lazarides  _et al._
(2015); Abanin  _et al._ (2016), prethermalization Kuwahara  _et al._ (2016);
Canovi  _et al._ (2016); Machado  _et al._ (2019), and many body quantum scars
Turner  _et al._ (2018); Michailidis  _et al._ (2020); Pizzi  _et al._ (2020);
Surace  _et al._ (2021).

Here as in Ref. Kshetrimayum  _et al._ , 2021 the MBL mechanism is considered.
The Floquet Hamiltonian with period $T$ is

$H_{\rm{F}}=\left\\{\begin{array}[]{ll}H_{\rm
MBL},&,\rm{when~{}}t\in[0,T/2),\\\ H_{\rm
f},&,\rm{when~{}}t\in[T/2,T).\end{array}\right.$ (5)

Here $H_{\rm MBL}$ is the same as in (2) and

$H_{f}=(2\pi/T-2\epsilon)\sum_{j}S^{x}_{j}.$ (6)

For $\epsilon=0$ its evolution operator, $e^{-iH_{\rm f}T/2}$, is a perfect
spin flip operator that is applied once per every period $T$. In general it is
an imperfect spin rotation with an angle deficit $\epsilon T$. An initial
state is again the Néel state in (3).

Figure 3 shows examples of stroboscopic time evolutions for ancillary
dimensions $d_{A}=2,5$ and corresponding disorder strengths $h=50,100$ that
turn out to be strong enough to stabilize an unambiguous time crystalline
stage. The same figure demonstrates robustness of these time crystals against
the spin-flip imperfection with $\epsilon=0.5$. On the other hand, when the
disorder is not strong enough then the time crystalline stage either cannot be
identified at all or disappears within the achievable evolution time, see the
left and right panels in Fig. 4, respectively. The disappearance is
accompanied by a worsening convergence with $D$ as might have been expected
when MBL is not effective enough.

##  V Conclusion

The neighborhood tensor update algorithm Dziarmaga (2021) was employed to
simulate unitary time evolution of the square-lattice antiferromagnetic
spin-$1/2$ Heisenberg model with discrete disorder. Starting from the Néel
state many body localized regimes were identified — for both the 2-level and
the 5-level disorder — where the system retains some memory of the initial
state. When promoted to a Floquet Hamiltonian with periodic (imperfect) spin
flips the simulations revealed extended time-crystalline stages for strong
enough disorder.

###### Acknowledgements.

I am indebted for stimulating comments from Marin Bukov, Piotr Czarnik, Marek
Rams, and Aritra Sinha. This research was supported in part by the National
Science Centre (NCN), Poland under project 2019/35/B/ST3/01028.

Figure 5:  Trotter gate.  In (a) a 2-site gate is applied to physical indices
of NN tensors $A$ and $B$ as in Fig. 1 (b). The gate is replaced by two
tensors, $G_{A}$ and $G_{B}$, contracted by an index with dimension $r$. In
(b) the tensor contraction $A\cdot G_{A}$ is QR-decomposed into $Q_{A}R_{A}$.
Similarly $B\cdot G_{B}=Q_{B}R_{B}$. Isometries $Q_{A,B}$ will remain fixed.
In (c) after SVD, $R_{A}R_{B}^{T}=U_{A}SU_{B}^{T}$, $S$ is truncated to $D$
leading singular values. In (d) matrices $M_{A}=U_{A}S^{1/2}$ and
$M_{B}^{T}=S^{1/2}U_{B}^{T}$ are made by absorbing a square root of truncated
$S$ symmetrically. In (e) at this point new iPEPS tensors could be obtained as
$A^{\prime}=Q_{A}\cdot M_{A}$ and $B^{\prime}=Q_{B}\cdot M_{B}$ ending the
story. This scheme was referred to as an SVD update (SVDU) in Ref. Dziarmaga,
2021. In NTU scheme matrices $M_{A,B}$ are further optimized in the
neighborhood tensor environment in Fig. 1 before being contracted back with
the fixed isometries $Q_{A,B}$ to make new iPEPS tensors $A^{\prime}$ and
$B^{\prime}$.  Figure 6:  Metric tensor.  Norm squared of the matrix product,
$||M_{A}M_{B}^{T}||^{2}$, is shown in (d). Here $g$ is a metric tensor
assembled in (c). The upper/lower pair of free indices in (c) corresponds to
the upper/lower pair of indices of $g$ in (d). Diagram (c) is assembled from 6
edge tensors – two examples of which are shown in (a) – and 2 double isometric
tensors in (b). Cost of all the contractions scales like $D^{8}$ and is fully
parallelizable.

##  Appendix A Algorithm

The core of NTU scheme Dziarmaga (2021) is outlined in figures 1, 5, 6, and 7.
Figure 5 shows application of a 2-site Trotter gate to a pair of NN iPEPS
tensors, $A$ and $B$. The rank-$r$ gate enlarges the bond dimension from $D$
to $Dr$ that has to be truncated back to $D$ in an optimal way. For the sake
of efficiency, a QR decomposition is introducing reduced matrices $R_{A}$ and
$R_{B}$ in place of full tensors Evenbly (2018). A SVD of their product,
$R_{A}R_{B}^{T}$, can be truncated to provide a preliminary truncation of the
bond dimension. The truncated reduced matrices, $M_{A}$ and $M_{B}$, are
further optimized iteratively in Figs. 6 and 7. After convergence they are
contracted again with the fixed isometries, $Q_{A}$ and $Q_{B}$, into new
iPEPS tensors $A^{\prime}$ and $B^{\prime}$. The numerical cost of the
procedures in Fig. 5 scales as $D^{5}$.

The further NTU optimization minimizes the difference between the LHS and the
RHS of the equation in Fig. 1(b). The contraction $A^{\prime}-B^{\prime}$ on
its RHS (the purple tensors) is the same contraction as in Fig. 5(e) and it
has to be understood as the diagram in Fig. 5(d). This way the RHS depends on
the product of the matrices to be optimized: $M_{A}M_{B}^{T}$. Therefore, the
norm squared of the difference between the LHS and the RHS can be written as

$F\left(M_{A}M_{B}^{T}\right)=\left[M_{A}M_{B}^{T}-R_{A}R_{B}^{T}\right]^{\dagger}~{}g~{}\left[M_{A}M_{B}^{T}-R_{A}R_{B}^{T}\right],$
(7)

where $g$ is the metric tensor constructed in Fig. 6. With metric $g$ fixed,
matrices $M_{A}$ and $M_{B}$ are optimized to make them the best approximation
to the untruncated/exact product $R_{A}R_{B}^{T}$. For a fixed $M_{B}$ the
cost function (7) becomes a quadratic form in $M_{A}$:

$F_{A}\left(M_{A}\right)=M_{A}^{\dagger}g_{A}M_{A}-M_{A}^{\dagger}J_{A}-J_{A}^{\dagger}M_{A}+F(0),$
(8)

where $g_{A}$ and $J_{A}$ depend on the fixed $M_{B}$, see Figs. 7(b) and (c).
The matrix is optimized as

$M_{A}={\rm pinv}\left(g_{A}\right)J_{A},$ (9)

where tolerance of the pseudo-inverse can be adjusted to minimize
$F_{A}\left[{\rm pinv}\left(g_{A}\right)J_{A}\right]$. Thanks to the exactness
of $g$ in NTU, the optimal tolerance is usually in the range
$10^{-12},...,10^{-8}$. This optimization of $M_{A}$ is followed by a similar
optimization of $M_{B}$. The optimizations are repeated in a loop,

$\rightarrow M_{A}\rightarrow M_{B}\rightarrow,$ (10)

until convergence of a relative NTU error:
$F\left(M_{A}M_{B}^{T}\right)/F(0)$. In general the converged error is non-
zero due to the necessary truncation of the bond dimension.

The converged error is used as a criterion to terminate the time evolution.
More specifically, we use

$\delta=dt^{-1}\sqrt{F\left(M_{A}M_{B}^{T}\right)/F(0)}$ (11)

to measure accuracy of the Trotterized evolution. It should be independent of
the time step for small enough $dt$. The square root makes $dt\cdot\delta$ a
rough estimator of a relative error of the wavefunction inflicted by the
Trotter gate and as such also of an error of its expectation values. In the
worst case scenario, where the gate errors accumulate in an additive way, a
typical relative error after time $t$ is proportional to $t\cdot\delta$, where
$\delta$ is averaged over the evolution time.

Figure 7:  NTU optimization loop.  In (a) matrices $M_{A},M_{B}$ are optimized
for their product, $M_{A}M_{B}^{T}$, to be the best approximation to the exact
product, $R_{A}R_{B}^{T}$. The error is measured with the metric in Fig. 6(d).
In (b) reduced metric tensor $g_{A}$ for matrix $M_{A}$. In (c) reduced source
term $J_{A}$ for matrix $M_{A}$. In (d) a product of converged matrices is
subject to a SVD, $M_{A}M_{B}=U_{A}SU_{B}^{T}$, after which new balanced
matrices, $M_{A}=U_{A}S^{1/2}$ and $M_{B}^{T}=S^{1/2}U_{B}^{T}$, are formed by
absorbing singular values $S$ in a symmetric way. However, iterative
optimization of the matrices is not symmetric. Before optimization with
respect to $M_{A}$ the matrices are “tilted” as $M_{A}=U_{A}S$ and
$M_{B}^{T}=U_{B}^{T}$ and vice versa Evenbly (2018).

Except for SVD of small matrices, $R_{A}R_{B}^{T}$ and $M_{A}M_{B}^{T}$, the
$\propto D^{8}$ operations – required to simulate time evolution in the NTU
scheme – are fully parallelizable. A potential bottleneck is calculation of
expectations values that requires the corner transfer matrix renormalization
group Orús (2014). However, the expectation values usually do not need to be
evaluated after every time step and they may not require the same precision as
stable evolution with the full update (FU) scheme.

## References

  * Verstraete  _et al._ (2008) F. Verstraete, V. Murg, and J. Cirac, Adv. Phys. 57, 143 (2008).
  * Orús (2014) R. Orús, Ann. Phys. (Amsterdam) 349, 117 (2014).
  * Fannes  _et al._ (1992) M. Fannes, B. Nachtergaele, and R. Werner, Comm. in Math. Phys. 144, 443 (1992).
  * Nishio  _et al._ (2004) Y. Nishio, N. Maeshima, A. Gendiar, and T. Nishino, arXiv:cond-mat/0401115 (2004).
  * Verstraete and Cirac (2004a) F. Verstraete and J. I. Cirac, arXiv:cond-mat/0407066 (2004a).
  * Vidal (2007) G. Vidal, Phys. Rev. Lett. 99, 220405 (2007).
  * Vidal (2008) G. Vidal, Phys. Rev. Lett. 101, 110501 (2008).
  * Evenbly and Vidal (2014a) G. Evenbly and G. Vidal, Phys. Rev. Lett. 112, 220502 (2014a).
  * Evenbly and Vidal (2014b) G. Evenbly and G. Vidal, Phys. Rev. B 89, 235113 (2014b).
  * Hastings (2007) M. B. Hastings, J. Stat. Mech. Theory Exp. 2007, P08024 (2007).
  * Schuch  _et al._ (2008) N. Schuch, M. M. Wolf, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 100, 030504 (2008).
  * Barthel (2017) T. Barthel, arXiv:1708.09349 (2017).
  * White (1992) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
  * White (1993) S. R. White, Phys. Rev. B 48, 10345 (1993).
  * Schollwöck (2005) U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).
  * Schöllwock (2011) U. Schöllwock, Ann. Phys. (Amsterdam) 326, 96 (2011).
  * Wolf  _et al._ (2008) M. M. Wolf, F. Verstraete, M. B. Hastings, and J. I. Cirac, Phys. Rev. Lett. 100, 070502 (2008).
  * Molnar  _et al._ (2015) A. Molnar, N. Schuch, F. Verstraete, and J. I. Cirac, Phys. Rev. B 91, 045138 (2015).
  * Ge and Eisert (2016) Y. Ge and J. Eisert, New J. Phys. 18, 083026 (2016).
  * Corboz  _et al._ (2010a) P. Corboz, G. Evenbly, F. Verstraete, and G. Vidal, Phys. Rev. A 81, 010303(R) (2010a).
  * Pineda  _et al._ (2010) C. Pineda, T. Barthel, and J. Eisert, Phys. Rev. A 81, 050303(R) (2010).
  * Corboz and Vidal (2009) P. Corboz and G. Vidal, Phys. Rev. B 80, 165129 (2009).
  * Barthel  _et al._ (2009) T. Barthel, C. Pineda, and J. Eisert, Phys. Rev. A 80, 042333 (2009).
  * Gu  _et al._ (2010) Z.-C. Gu, F. Verstraete, and X.-G. Wen, arXiv:1004.2563 (2010).
  * Kraus  _et al._ (2010) C. V. Kraus, N. Schuch, F. Verstraete, and J. I. Cirac, Phys. Rev. A 81, 052338 (2010).
  * Corboz  _et al._ (2010b) P. Corboz, R. Orús, B. Bauer, and G. Vidal, Phys. Rev. B 81, 165104 (2010b).
  * Corboz  _et al._ (2011) P. Corboz, S. R. White, G. Vidal, and M. Troyer, Phys. Rev. B 84, 041108(R) (2011).
  * Verstraete and Cirac (2004b) F. Verstraete and J. I. Cirac, cond-mat/0407066 (2004b).
  * Murg  _et al._ (2007) V. Murg, F. Verstraete, and J. I. Cirac, Phys. Rev. A 75, 033605 (2007).
  * Jordan  _et al._ (2008) J. Jordan, R. Orús, G. Vidal, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 101, 250602 (2008).
  * Jiang  _et al._ (2008) H. C. Jiang, Z. Y. Weng, and T. Xiang, Phys. Rev. Lett. 101, 090603 (2008).
  * Gu  _et al._ (2008) Z.-C. Gu, M. Levin, and X.-G. Wen, Phys. Rev. B 78, 205116 (2008).
  * Orús and Vidal (2009) R. Orús and G. Vidal, Phys. Rev. B 80, 094403 (2009).
  * Matsuda  _et al._ (2013) Y. H. Matsuda, N. Abe, S. Takeyama, H. Kageyama, P. Corboz, A. Honecker, S. R. Manmana, G. R. Foltin, K. P. Schmidt, and F. Mila, Phys. Rev. Lett. 111, 137204 (2013).
  * Corboz and Mila (2014) P. Corboz and F. Mila, Phys. Rev. Lett. 112, 147203 (2014).
  * Zheng  _et al._ (2017) B.-X. Zheng, C.-M. Chung, P. Corboz, G. Ehlers, M.-P. Qin, R. M. Noack, H. Shi, S. R. White, S. Zhang, and G. K.-L. Chan, Science 358, 1155 (2017).
  * Liao  _et al._ (2017) H. J. Liao, Z. Y. Xie, J. Chen, Z. Y. Liu, H. D. Xie, R. Z. Huang, B. Normand, and T. Xiang, Phys. Rev. Lett. 118, 137202 (2017).
  * Phien  _et al._ (2015) H. N. Phien, J. A. Bengua, H. D. Tuan, P. Corboz, and R. Orús, Phys. Rev. B 92, 035142 (2015).
  * Corboz (2016a) P. Corboz, Phys. Rev. B 94, 035133 (2016a).
  * Vanderstraeten  _et al._ (2016) L. Vanderstraeten, J. Haegeman, P. Corboz, and F. Verstraete, Phys. Rev. B 94, 155123 (2016).
  * Fishman  _et al._ (2018) M. T. Fishman, L. Vanderstraeten, V. Zauner-Stauber, J. Haegeman, and F. Verstraete, Phys. Rev. B 98, 235148 (2018).
  * Xie  _et al._ (2017) Z. Y. Xie, H. J. Liao, R. Z. Huang, H. D. Xie, J. Chen, Z. Y. Liu, and T. Xiang, Phys. Rev. B 96, 045128 (2017).
  * Corboz (2016b) P. Corboz, Phys. Rev. B 93, 045116 (2016b).
  * Corboz  _et al._ (2018) P. Corboz, P. Czarnik, G. Kapteijns, and L. Tagliacozzo, Phys. Rev. X 8, 031031 (2018).
  * Rader and Läuchli (2018) M. Rader and A. M. Läuchli, Phys. Rev. X 8, 031030 (2018).
  * Rams  _et al._ (2018) M. M. Rams, P. Czarnik, and L. Cincio, Phys. Rev. X 8, 041033 (2018).
  * Czarnik  _et al._ (2012) P. Czarnik, L. Cincio, and J. Dziarmaga, Phys. Rev. B 86, 245101 (2012).
  * Czarnik and Dziarmaga (2014) P. Czarnik and J. Dziarmaga, Phys. Rev. B 90, 035144 (2014).
  * Czarnik and Dziarmaga (2015a) P. Czarnik and J. Dziarmaga, Phys. Rev. B 92, 035120 (2015a).
  * Czarnik  _et al._ (2016a) P. Czarnik, J. Dziarmaga, and A. M. Oleś, Phys. Rev. B 93, 184410 (2016a).
  * Czarnik and Dziarmaga (2015b) P. Czarnik and J. Dziarmaga, Phys. Rev. B 92, 035152 (2015b).
  * Czarnik  _et al._ (2016b) P. Czarnik, M. M. Rams, and J. Dziarmaga, Phys. Rev. B 94, 235142 (2016b).
  * Czarnik  _et al._ (2017) P. Czarnik, J. Dziarmaga, and A. M. Oleś, Phys. Rev. B 96, 014420 (2017).
  * Dai  _et al._ (2017) Y.-W. Dai, Q.-Q. Shi, S. Y. Cho, M. T. Batchelor, and H.-Q. Zhou, Phys. Rev. B 95, 214409 (2017).
  * Czarnik  _et al._ (2019a) P. Czarnik, J. Dziarmaga, and P. Corboz, Phys. Rev. B 99, 035115 (2019a).
  * Czarnik and Corboz (2019) P. Czarnik and P. Corboz, Phys. Rev. B 99, 245107 (2019).
  * Kshetrimayum  _et al._ (2019) A. Kshetrimayum, M. Rizzi, J. Eisert, and R. Orús, Phys. Rev. Lett. 122, 070502 (2019).
  * Czarnik  _et al._ (2019b) P. Czarnik, A. Francuz, and J. Dziarmaga, Phys. Rev. B 100, 165147 (2019b).
  * Jiménez  _et al._ (2021) J. L. Jiménez, S. P. G. Crone, E. Fogh, M. E. Zayed, R. Lortz, E. Pomjakushina, K. Conder, A. M. Läuchli, L. Weber, S. Wessel, A. Honecker, B. Normand, C. Rüegg, P. Corboz, H. M. Rønnow, and F. Mila, Nature 592, 370 (2021).
  * Jiménez  _et al._ (2020) J. L. Jiménez, S. P. G. Crone, E. Fogh, M. E. Zayed, R. Lortz, E. Pomjakushina, K. Conder, A. M. Läuchli, L. Weber, S. Wessel, A. Honecker, B. Normand, C. Rüegg, P. Corboz, H. M. Rønnow, and F. Mila, arXiv:2009.14492 (2020).
  * Poilblanc  _et al._ (2020) D. Poilblanc, M. Mambrini, and F. Alet, arXiv:2010.07828 (2020).
  * Czarnik  _et al._ (2021) P. Czarnik, M. M. Rams, P. Corboz, and J. Dziarmaga, Phys. Rev. B 103, 075113 (2021).
  * Poilblanc  _et al._ (2021) D. Poilblanc, M. Mambrini, and F. Alet, SciPost Phys. 10, 19 (2021).
  * Kshetrimayum  _et al._ (2017) A. Kshetrimayum, H. Weimer, and R. Orús, Nat. Commun. 8, 1291 (2017).
  * Mc Keever and Szymańska (2021) C. Mc Keever and M. H. Szymańska, Phys. Rev. X 11, 021035 (2021).
  * Vanderstraeten  _et al._ (2015) L. Vanderstraeten, M. Mariën, F. Verstraete, and J. Haegeman, Phys. Rev. B 92, 201111(R) (2015).
  * Ponsioen and Corboz (2020) B. Ponsioen and P. Corboz, Phys. Rev. B 101, 195109 (2020).
  * Hubig and Cirac (2019) C. Hubig and J. I. Cirac, SciPost Phys. 6, 31 (2019).
  * Hubig  _et al._ (2020) C. Hubig, A. Bohrdt, M. Knap, F. Grusdt, and J. I. Cirac, SciPost Phys. 8, 21 (2020).
  * Abendschein and Capponi (2008) A. Abendschein and S. Capponi, Phys. Rev. Lett. 101, 227201 (2008).
  * Kshetrimayum  _et al._ (2020) A. Kshetrimayum, M. Goihl, and J. Eisert, Phys. Rev. B 102, 235132 (2020).
  * Kshetrimayum  _et al._ (2021) A. Kshetrimayum, M. Goihl, D. M. Kennes, and J. Eisert, Phys. Rev. B 103, 224205 (2021).
  * Dziarmaga (2021) J. Dziarmaga, Phys. Rev. B 104, 094411 (2021).
  * Schmitt  _et al._ (2021) M. Schmitt, M. M. Rams, J. Dziarmaga, M. Heyl, and W. H. Zurek, “Quantum phase transition dynamics in the two-dimensional transverse-field ising model,” (2021), arXiv:2106.09046 [cond-mat.str-el] .
  * Kaneko and Danshita (2021) R. Kaneko and I. Danshita, “Tensor-network study of correlation-spreading dynamics in the two-dimensional bose-hubbard model,” (2021), arXiv:2108.11051 [cond-mat.quant-gas] .
  * Cincio and Vidal (2013) L. Cincio and G. Vidal, Phys. Rev. Lett. 110, 067208 (2013).
  * Bruognolo  _et al._ (2017) B. Bruognolo, Z. Zhu, S. R. White, and E. M. Stoudenmire, arXiv:1705.05578 (2017).
  * Chen  _et al._ (2018a) B.-B. Chen, L. Chen, Z. Chen, W. Li, and A. Weichselbaum, Phys. Rev. X 8, 031082 (2018a).
  * Chen  _et al._ (2019) L. Chen, D.-W. Qu, H. Li, B.-B. Chen, S.-S. Gong, J. von Delft, A. Weichselbaum, and W. Li, Phys. Rev. B 99, 140404(R) (2019).
  * Li  _et al._ (2019) H. Li, B.-B. Chen, Z. Chen, J. von Delft, A. Weichselbaum, and W. Li, Phys. Rev. B 100, 045110 (2019).
  * Chen  _et al._ (2020) B.-B. Chen, C. Chen, Z. Chen, J. Cui, Y. Zhai, A. Weichselbaum, J. von Delft, Z. Y. Meng, and W. Li, arXiv:2008.02179 (2020).
  * Li  _et al._ (2011) W. Li, S.-J. Ran, S.-S. Gong, Y. Zhao, B. Xi, F. Ye, and G. Su, Phys. Rev. Lett. 106, 127202 (2011).
  * Xie  _et al._ (2012) Z. Y. Xie, J. Chen, M. P. Qin, J. W. Zhu, L. P. Yang, and T. Xiang, Phys. Rev. B 86, 045139 (2012).
  * Ran  _et al._ (2012) S.-J. Ran, W. Li, B. Xi, Z. Zhang, and G. Su, Phys. Rev. B 86, 134429 (2012).
  * Ran  _et al._ (2013) S.-J. Ran, B. Xi, T. Liu, and G. Su, Phys. Rev. B 88, 064407 (2013).
  * Ran  _et al._ (2018) S.-J. Ran, W. Li, S.-S. Gong, A. Weichselbaum, J. von Delft, and G. Su, Phys. Rev. B 97, 075146 (2018).
  * Peng  _et al._ (2017) C. Peng, S.-J. Ran, T. Liu, X. Chen, and G. Su, Phys. Rev. B 95, 075140 (2017).
  * Chen  _et al._ (2018b) X. Chen, S.-J. Ran, T. Liu, C. Peng, Y.-Z. Huang, and G. Su, Science Bulletin 63, 1545 (2018b).
  * Ran  _et al._ (2019) S.-J. Ran, B. Xi, C. Peng, G. Su, and M. Lewenstein, Phys. Rev. B 99, 205132 (2019).
  * Vanhecke  _et al._ (2021) B. Vanhecke, D. Devoogdt, F. Verstraete, and L. Vanderstraeten, “Simulating thermal density operators with cluster expansions and tensor networks,” (2021), arXiv:2112.01507 [cond-mat.str-el] .
  * Wang and Verstraete (2011) L. Wang and F. Verstraete, (2011), arXiv:1110.4362 [cond-mat.str-el] .
  * Lubasch  _et al._ (2014a) M. Lubasch, J. I. Cirac, and M.-C. Bañuls, New Journal of Physics 16, 033014 (2014a).
  * Lubasch  _et al._ (2014b) M. Lubasch, J. I. Cirac, and M.-C. Bañuls, Phys. Rev. B 90, 064425 (2014b).
  * Basko  _et al._ (2006) D. Basko, I. Aleiner, and B. Altshuler, Annals of Physics 321, 1126 (2006).
  * Abanin  _et al._ (2019) D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev. Mod. Phys. 91, 021001 (2019).
  * Pal and Huse (2010) A. Pal and D. A. Huse, Phys. Rev. B 82, 174411 (2010).
  * Vosk  _et al._ (2015) R. Vosk, D. A. Huse, and E. Altman, Phys. Rev. X 5, 031032 (2015).
  * Goihl  _et al._ (2018) M. Goihl, M. Gluza, C. Krumnow, and J. Eisert, Phys. Rev. B 97, 134202 (2018).
  * Bordia  _et al._ (2017) P. Bordia, H. Lüschen, S. Scherg, S. Gopalakrishnan, M. Knap, U. Schneider, and I. Bloch, Phys. Rev. X 7, 041047 (2017).
  * De Roeck and Huveneers (2017) W. De Roeck and F. m. c. Huveneers, Phys. Rev. B 95, 155129 (2017).
  * Geißler and Pupillo (2020) A. Geißler and G. Pupillo, Phys. Rev. Research 2, 042037 (2020).
  * Wahl  _et al._ (2019) T. B. Wahl, A. Pal, and S. H. Simon, Nature Physics 15, 164 (2019).
  * Kennes (2018) D. M. Kennes, “Many-body localization in two dimensions from projected entangled-pair states,” (2018), arXiv:1811.04126 [cond-mat.dis-nn] .
  * De Tomasi  _et al._ (2019) G. De Tomasi, F. Pollmann, and M. Heyl, Phys. Rev. B 99, 241114 (2019).
  * Andraschko  _et al._ (2014) F. Andraschko, T. Enss, and J. Sirker, Phys. Rev. Lett. 113, 217201 (2014).
  * Paredes  _et al._ (2005) B. Paredes, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 95, 140501 (2005).
  * Wilczek (2012) F. Wilczek, Phys. Rev. Lett. 109, 160401 (2012).
  * Watanabe and Oshikawa (2015) H. Watanabe and M. Oshikawa, Phys. Rev. Lett. 114, 251603 (2015).
  * Kozin and Kyriienko (2019) V. K. Kozin and O. Kyriienko, Phys. Rev. Lett. 123, 210602 (2019).
  * Sacha (2015) K. Sacha, Phys. Rev. A 91, 033617 (2015).
  * Khemani  _et al._ (2016) V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. 116, 250401 (2016).
  * Else  _et al._ (2016) D. V. Else, B. Bauer, and C. Nayak, Phys. Rev. Lett. 117, 090402 (2016).
  * Yao  _et al._ (2017) N. Y. Yao, A. C. Potter, I.-D. Potirniche, and A. Vishwanath, Phys. Rev. Lett. 118, 030401 (2017).
  * Zhang  _et al._ (2017) J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I.-D. Potirniche, A. C. Potter, A. Vishwanath, N. Y. Yao, and C. Monroe, Nature 543, 217 (2017).
  * Choi  _et al._ (2017) S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, C. von Keyserlingk, N. Y. Yao, E. Demler, and M. D. Lukin, Nature 543, 221 (2017).
  * Rovny  _et al._ (2018a) J. Rovny, R. L. Blum, and S. E. Barrett, Phys. Rev. Lett. 120, 180603 (2018a).
  * Pal  _et al._ (2018) S. Pal, N. Nishad, T. S. Mahesh, and G. J. Sreejith, Phys. Rev. Lett. 120, 180602 (2018).
  * Rovny  _et al._ (2018b) J. Rovny, R. L. Blum, and S. E. Barrett, Phys. Rev. B 97, 184301 (2018b).
  * Smits  _et al._ (2018) J. Smits, L. Liao, H. T. C. Stoof, and P. van der Straten, Phys. Rev. Lett. 121, 185301 (2018).
  * Liao  _et al._ (2019) L. Liao, J. Smits, P. van der Straten, and H. T. C. Stoof, Phys. Rev. A 99, 013625 (2019).
  * Smits  _et al._ (2020) J. Smits, H. T. C. Stoof, and P. van der Straten, New Journal of Physics 22, 105001 (2020).
  * Sacha and Zakrzewski (2017) K. Sacha and J. Zakrzewski, Reports on Progress in Physics 81, 016401 (2017).
  * Khemani  _et al._ (2019) V. Khemani, R. Moessner, and S. L. Sondhi, “A brief history of time crystals,” (2019), arXiv:1910.10745 [cond-mat.str-el] .
  * Else  _et al._ (2020) D. V. Else, C. Monroe, C. Nayak, and N. Y. Yao, Annual Review of Condensed Matter Physics 11, 467 (2020), https://doi.org/10.1146/annurev-conmatphys-031119-050658 .
  * Ponte  _et al._ (2015) P. Ponte, Z. Papić, F. m. c. Huveneers, and D. A. Abanin, Phys. Rev. Lett. 114, 140401 (2015).
  * Lazarides  _et al._ (2015) A. Lazarides, A. Das, and R. Moessner, Phys. Rev. Lett. 115, 030402 (2015).
  * Abanin  _et al._ (2016) D. A. Abanin, W. De Roeck, and F. Huveneers, Annals of Physics 372, 1 (2016).
  * Kuwahara  _et al._ (2016) T. Kuwahara, T. Mori, and K. Saito, Annals of Physics 367, 96 (2016).
  * Canovi  _et al._ (2016) E. Canovi, M. Kollar, and M. Eckstein, Phys. Rev. E 93, 012130 (2016).
  * Machado  _et al._ (2019) F. Machado, G. D. Kahanamoku-Meyer, D. V. Else, C. Nayak, and N. Y. Yao, Phys. Rev. Research 1, 033202 (2019).
  * Turner  _et al._ (2018) C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Nature Physics 14, 745 (2018).
  * Michailidis  _et al._ (2020) A. A. Michailidis, C. J. Turner, Z. Papić, D. A. Abanin, and M. Serbyn, Phys. Rev. X 10, 011055 (2020).
  * Pizzi  _et al._ (2020) A. Pizzi, D. Malz, G. De Tomasi, J. Knolle, and A. Nunnenkamp, Phys. Rev. B 102, 214207 (2020).
  * Surace  _et al._ (2021) F. M. Surace, M. Votto, E. G. Lazo, A. Silva, M. Dalmonte, and G. Giudici, Phys. Rev. B 103, 104302 (2021).
  * Evenbly (2018) G. Evenbly, Phys. Rev. B 98, 085155 (2018).


# Is there a classical model of Wigner’s friend?

Anthony Sudbery1
Department of Mathematics, University of York,
Heslington, York, England YO10 5DD
1<EMAIL_ADDRESS>
(31 October 2021)

###### Abstract

“Wigner’s friend” refers to a quantum process of which different observers,
following the rules of quantum mechanics, give contradictory descriptions.
Lostaglio and Bowles have recently claimed to describe a classical system
showing the same effect. It is argued that this claim is not justified. It
fails to take account of the different meanings of probability in quantum and
classical mechanics.

In a famous discussion of the foundations of quantum mechanics [3], Wigner
described a situation in which two observers use quantum mechanics to describe
the same physical state of affairs, but their descriptions are not compatible.
The first observer — Wigner’s friend, who I will call Frieda — takes a qubit
in the state $\tfrac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$ and measures it in the
basis $\\{|0\rangle,|1\rangle\\}$, observing the result $0$ or $1$; let us
suppose it is $0$. Following the directions in her quantum mechanics textbook,
she concludes that she and the qubit are together in a state $|00\rangle$, in
an obvious notation. Wigner, however, observing Frieda and the lab from
outside, sees that the effect of the measurement is to entangle the qubit with
Frieda, and concludes that the resulting state of them both is
$|+\rangle=\tfrac{1}{\sqrt{2}}(|00\rangle+|11\rangle)$.

These descriptions are not compatible; they have implications that are clearly
contradictory. Wigner can confirm his description of the state of the lab by
measuring it in the basis
$|\pm\rangle=\tfrac{1}{\sqrt{2}}(|00\rangle\pm|11\rangle)$; it follows from
his description that the result of the measurement will certainly be $+$.
Frieda, however, describes the state as
$|00\rangle=\tfrac{1}{\sqrt{2}}(|+\rangle+|-\rangle)$, from which it follows
that the result of Wigner’s measurement is equally likely to be $|+\rangle$ or
$|-\rangle$. This statement of probabilities can only be checked by
considering an ensemble of systems, i.e. by repeating the experiment in an
identical manner. In each repetition the state $|00\rangle$ can lead to either
result, and in a run of $N$ repetitions of the experiment all $2^{N}$
sequences of results are equally likely; if $N$ is large enough the result
“$-$” is practically certain to occur. It is a consequence of Wigner’s state
$|+\rangle$, however, that the result “$-$” will never occur.

This incompatibility can be described formally by considering the logical
conjunction of Frieda’s state and Wigner’s state. This conjunction, the
intersection of the two state vectors, is the zero vector, which represents
impossibility.

Lostaglio and Bowles [1] claim that this situation can be replicated in a
deterministic classical model. In this model the system consists of a set of
four bits $X,X^{\prime},Y,Y^{\prime}$ which always have definite values $0$ or
$1$, and states of the system are identified with probability distributions on
the sample space $\\{0,1\\}^{4}$, which can be written as a vector in
$\mathbb{R}^{4}$ which is a convex combination of the sixteen extreme points
$(\epsilon_{1})_{X}(\epsilon_{2})_{X^{\prime}}(\epsilon_{3})_{Y}(\epsilon_{4})_{Y^{\prime}}$
where $\epsilon_{i}\in\\{0,1\\}$. Such a state is a state of knowledge of an
agent; the probabilities (coefficients in the convex combination) measure the
strength of the agent’s belief in the values of the bits. Lostaglio and Bowles
describe a process which leads Wigner and Frieda to assign the following
states to the bits. Frieda’s state is either

$\tfrac{1}{4}\big{(}0_{X}0_{X^{\prime}}0_{Y}0_{Y^{\prime}}+0_{X}0_{X^{\prime}}0_{Y}1_{Y^{\prime}}+0_{X}1_{X^{\prime}}0_{Y}0_{Y^{\prime}}+0_{X}1_{X^{\prime}}0_{Y}1_{Y^{\prime}}\big{)}$
(1)

or

$\tfrac{1}{4}\big{(}1_{X}0_{X^{\prime}}1_{Y}0_{Y^{\prime}}+1_{X}0_{X^{\prime}}1_{Y}1_{Y^{\prime}}+1_{X}1_{X^{\prime}}1_{Y}0_{Y^{\prime}}+1_{X}1_{X^{\prime}}1_{Y}1_{Y^{\prime}}\big{)},$
(2)

depending on the result of her measurement. Wigner’s state is

$\tfrac{1}{4}\big{(}0_{X}0_{X^{\prime}}0_{Y}0_{Y^{\prime}}+0_{X}1_{X^{\prime}}0_{Y}1_{Y^{\prime}}+1_{X}0_{X^{\prime}}1_{Y}1_{Y^{\prime}}+1_{X}1_{X^{\prime}}1_{Y}0_{Y^{\prime}}\big{)}.$
(3)

Unlike the state assignments in the quantum model, there is no contradiction
in the above descriptions. The descriptions are different, for the simple
reason that Wigner and Frieda have different information, and different areas
of ignorance; but in both cases (of Frieda’s result) they are compatible, not
contradictory. Logically, disjunctions “P or Q or R or S” and “P or S or T or
U” do not contradict each other, and this remains true if the disjunctions
have epistemic probabilities attached to them. If Wigner and Frieda were to
compare their knowledge, they would both update their states — to
$\tfrac{1}{2}\big{(}0_{X}0_{X^{\prime}}0_{Y}0_{Y^{\prime}}+0_{X}1_{X^{\prime}}0_{Y}1_{Y^{\prime}}\big{)}$
in case (1), to
$\tfrac{1}{2}\big{(}1_{X}0_{X^{\prime}}1_{Y}1_{Y^{\prime}}+1_{X}1_{X^{\prime}}1_{Y}0_{Y^{\prime}}\big{)}$
in case (2).

Lostaglio and Bowles now represent a measurement on the 4-bit system by a
deterministic process which takes each of the above states to a definite state
from which the values of the two bits $X$ and $Y$ (Frieda and her qubit) can
be predicted. The result is that Wigner predicts the values $(0,0)$ with
probability 1, while Frieda, in both cases, predicts that the values will be
either $(0,0)$ or $(1,0)$ with probability $1/2$ each. Again, these
predictions are different, but they are not contradictory. The logical
conjunction of the two confirms Wigner’s prediction of $(0,0)$.

To see why the situation is different in the two models, classical and
quantum, consider, as we did in the quantum case, what would happen if the
experiment was repeated a large number of times. Because the four bits have
definite values which will be the same in each repetition of the experiment,
and the process is deterministic, the result will be the same in every
repetition. Frieda’s prediction is that in a run of $N$ experiments there are
just two possible results: either the result is $(0,0)$ every time, or it is
$(1,0)$ every time. Her assignment of probability $1/2$ to the result $(1,0)$
does not imply that that result is bound to occur, contradicting Wigner. In
the quantum case, by contrast, the process of measurement is stochastic, the
repetitions are independent, and there are $2^{N}$ possible results; if $N$ is
large enough the result $|-\rangle$, because it has probability $1/2$, is
practically certain to occur in one of the repetitions.

The difference arises because the probabilities of quantum mechanics refer to
future measurements, and ensembles of such measurements, and these
probabilities concern what actually happens stochastically in the real world.
The probabilities in Lostaglio and Bowles’s model, on the other hand, are
measures of present knowledge, as in thermodynamics, and the measurement is
deterministic. As shown above, in the quantum situation the predictions
obtained from Wigner’s state and Frieda’s state do contradict each other, and
if they were to compare their states, they would have no way of resolving
their disagreement. There is a genuine paradox. This is true whether the
nature of the states is epistemic or ontic. The conclusion must be that the
theory which generates these predictions — universally valid quantum mechanics
with the projection postulate applied to all measurements — is untenable. It
seems unlikely that any model which might be called “classical” could exhibit
such contradictions.

## References

  * [1] Matteo Lostaglio and Joseph Bowles.  The original Wigner’s friend paradox within a realist toy model.  Proc. Roy. Soc. A., 477:20210273, 2021.  arXiv:2102.11032.
  * [2] J. A. Wheeler and W. H. Zurek, editors.  Quantum Theory and Measurement.  Princeton University Press, 1983.
  * [3] E. P. Wigner.  Remarks on the mind-body question.  In I. J. Good, editor, The Scientist Speculates, pages 284–302. Heinemann, 1962.  reprinted in [2] pp. 168-181.


#  Coulomb phase corrections to the transverse analyzing power
$A_{\text{N}}(t)$ in high energy forward proton-proton scattering 111This
manuscript has been authored by employees of Brookhaven Science Associates,
LLC under Contract No. DE-SC0012704 with the U.S. Department of Energy

A. A. Poblaguev<EMAIL_ADDRESS>Brookhaven National Laboratory, Upton, New
York 11973, USA

(November 2, 2021)

###### Abstract

Study of polarized proton-proton elastic scattering in the Coulomb-nuclear
interference region allows one to measure the forward hadronic single spin-
flip amplitude including its phase. However, in a precision experimental data
analysis, a phase shift correction $\delta_{C}$ due to the long distance
Coulomb interaction should be taken into account. For unpolarized scattering,
$\delta_{C}$ is commonly considered as well established. Here, we evaluate the
Coulomb phase shifts for the forward elastic proton-proton single spin-flip
electromagnetic and hadronic amplitudes. Only a small discrepancy between the
spin-flip and non-flip phases was found which can be neglected in the high
energy forward elastic $\mathit{pp}$ studies involving transverse spin.
Nonetheless, the effective alteration of the hadronic spin-flip amplitude by
the long distance electromagnetic corrections can be essential for
interpretation of the experimental results.

###### keywords:

Coulomb nuclear interference; Analyzing power; Coulomb phase shift; Absorptive
corrections;

††journal: Physics Letters B

##  1 Introduction

It is well known that a hadronic amplitude can be effectively altered by the
accompanying electromagnetic interaction. In the Coulomb-nuclear interference
(CNI) region, the elastic proton-proton amplitude can be approximated by the
following sum of hadronic $\phi(s,t)$ and electromagnetic, $t_{c}/t$, ones [1]

$\phi_{pp}^{\text{CNI}}(s,t)=\text{Im}\,\phi(s,0)\left[(i\\!+\\!\rho)e^{Bt/2}\\!+\\!\frac{t_{c}}{t}e^{i\delta_{C}+\widetilde{B}t/2}\right]\\!.$
(1)

Here, $\delta_{C}(s,t)$ stands for difference between the electromagnetic and
hadronic phases induced by the long-distance Coulomb interaction. Following
QED calculations and using the optical theorem, one finds
$t_{c}\\!=\\!-8\pi\alpha/\sigma_{\text{tot}}(s)$ [2], where $\alpha$ is the
fine structure constant and $\sigma_{\text{tot}}$ is the total $pp$ cross
section. Generally, the discussed amplitudes are functions of total energy
squared $s$ and momentum transfer squared $t$.

In this paper, numerical estimates will be done for a 100 GeV proton beam
(typical for the Relativistic Heavy Ion Collider or the future Electron Ion
Collider) scattering on a fixed proton target. Therefore, $\rho\\!=\\!-0.079$
[3], $\sigma_{\text{tot}}\\!=\\!39.2\,\text{mb}$ [3],
$t_{c}\\!=\\!-1.86\\!\times\\!10^{-3}\,\text{GeV}^{2}$, and
$B\\!=\\!11.2\,\text{GeV}^{-2}$ [4, 5]. For the sake of simplicity, we do not
distinguish between the hadronic $B$ and electromagnetic $\widetilde{B}$
slopes in expressions for the elastic differential cross section

$\frac{d\sigma^{el}}{dt}\propto\left[\left(\frac{t_{c}}{t}\right)^{2}-2\left(\rho+\delta_{C}\right)\frac{t_{c}}{t}+1+\rho^{2}\right]e^{Bt},$
(2)

and the analyzing power $A_{\text{N}}(t)$ [Eq. (8)]. We also do not consider
small corrections [6] due to Dirac and Pauli form factors and due to the
absorption [7].

For the unpolarized scattering, a theoretical understanding of $\delta_{C}(t)$
was developed in many works, particularly in [8, 9, 10, 11].

In Ref. [10], the Coulomb phase was evaluated as

$\delta_{C}=-\alpha\left[\ln{\frac{(B\\!+\\!\widetilde{B})|t|}{2}}+\gamma-\frac{\widetilde{B}t}{2}\ln{\frac{\widetilde{B}|t|}{2}}+\frac{\widetilde{B}t}{4}\right],$
(3)

where $\gamma\\!=\\!0.5772$ is Euler’s constant and
$\widetilde{B}\\!=\\!8/\Lambda^{2}\\!=\\!11.2\,\text{GeV}^{-2}$ was derived
from the $ep$ electromagnetic form factor in dipole form, ${\cal
F}_{ep}(t)\\!=\\!(1\\!-\\!t/\Lambda^{2})^{-2}$,
$\Lambda^{2}\\!=\\!0.71\,\text{GeV}^{2}$ [12]. Approach (3) is commonly used
in experimental data analysis.

Neglecting the slopes $B$ and $\widetilde{B}$ dependence on $t$ and omitting
the ${\cal O}(\alpha^{2})$ terms, an analytical solution for $\delta_{C}$ was
found in Ref. [11]

$\displaystyle\delta_{C}=\alpha\times\bigg{\\{}$
$\displaystyle\ln{\frac{\widetilde{B}^{2}}{B^{2}}}+Ei\left(\frac{2B^{2}/\widetilde{B}}{B+\widetilde{B}}w\right)-Ei(w)$
$\displaystyle+e^{2w}\left[2E_{1}(2w)-E_{1}(w)\right]\bigg{\\}},$ (4)

where $\mathit{Ei}(x)$ and $E_{1}(x)$ are the exponential integral functions,
$w\\!=\\!\left|\widetilde{B}t\right|/4$, and

$\widetilde{B}=\frac{2}{3}r_{E}^{2}\\!=\\!12.1\,\text{GeV}^{-2}$ (5)

was expressed via rms charge radius of a proton,
$r_{E}\\!=\\!0.841\,\text{fm}$ [13]. In terms of $B$ and $\widetilde{B}$, Eqs.
(3) and (4) are in agreement for low $t$.

It was shown in Ref. [14] that the Coulomb phase should be independent of the
helicity structure of the experimentally measured scattering amplitudes though
subsequent determination of the pure hadronic amplitudes may involve
order-$\alpha$ corrections resulting from the spin of the particles involved.

In this paper, we follow formulas derived in Refs. [11, 7] to evaluate Coulomb
phases in the single spin-flip scattering. The results obtained will be used
to evaluate the consequent alterations of the expression for forward elastic
${pp}$ analyzing power $A_{\text{N}}(t)$ suggested in Ref. [15].

##  2 High energy forward elastic $\boldsymbol{pp}$ analyzing power

For a CNI elastic scattering of a vertically polarized proton beam
$p^{\uparrow}{p}$, the analyzing power $A_{N}$ is defined by interference of
the non-flip (nf) and spin-flip (sf) amplitudes [2, 14, 15]

$A_{N}=\frac{2\,\text{Im}\left[\widetilde{\phi}_{\textit{sf}}\phi_{\textit{nf}}^{*}+\phi_{\textit{sf}}\widetilde{\phi}_{\textit{nf}}^{*}+\phi_{\textit{sf}}\phi_{\textit{nf}}^{*}\right]}{\left|\phi_{\textit{nf}}+\widetilde{\phi}_{\textit{nf}}\right|^{2}}.$
(6)

Here, hadronic $\phi$ and electromagnetic $\widetilde{\phi}$ parts of an
amplitude are discriminated by tilde symbol.

For $t\\!\to\\!0$, the sf amplitudes can be simply related to the nf ones:

$\widetilde{\phi}_{\textit{sf}}/\widetilde{\phi}_{\textit{nf}}=\frac{\kappa_{p}}{2}\,\sqrt{-t}/m_{p},\quad\phi_{\textit{sf}}/\phi_{\textit{nf}}=\frac{r_{5}}{i+\rho}\,\sqrt{-t}/m_{p},$
(7)

where $\kappa_{p}\\!=\\!\mu_{p}-1\\!=\\!1.973$ is anomalous magnetic moment of
a proton, $|r_{5}|\\!=\\!|R_{5}\\!+\\!iI_{5}|\\!\sim\\!0.02$ [16] is hadronic
spin-flip parameter [15], and $m_{p}$ is a proton mass. Consequently,

$\displaystyle A_{\text{N}}$ $\displaystyle(t)=\frac{\sqrt{-t}}{m_{p}}\times$
$\displaystyle\frac{\left[\kappa_{p}(1\\!-\\!\delta_{C}^{em}\rho)\\!-\\!2(I_{5}\\!-\\!\delta_{C}^{h}R_{5})\right]{t_{c}/t}\\!-\\!2(R_{5}\\!+\rho{I_{5}})}{\left({t_{c}/t}\right)^{2}-2\left(\rho+\delta_{C}\right){t_{c}/t}+1+\rho^{2}}\\!,$
(8)

where $\delta_{C}$ is non-flip phase (4) while $\delta_{C}^{em}$ and
$\delta_{C}^{h}$ are spin-flip phase shifts in the
$\widetilde{\phi}_{\textit{sf}}\phi_{\textit{nf}}^{*}$ and
$\phi_{\textit{sf}}\widetilde{\phi}_{\textit{nf}}^{*}$ interference terms,
respectively.

All recent experimental studies [16, 17, 18, 19] of the forward elastic
proton-proton $A_{\text{N}}(t)$ had been done assuming
$\delta_{C}^{em}\\!=\\!\delta_{C}^{h}\\!=\\!\delta_{C}$.

##  3 Spin-flip Coulomb phases in elastic
$\boldsymbol{{p^{\uparrow}}\\!\\!\\!{p}}$ analyzing power

Figure 1:  Three types of the elastic ${pp}$ scattering: (C) electromagnetic
including multi-photon exchange, (N) bare hadronic, and (NC) combined hadronic
and electromagnetic.

To evaluate spin-flip phases, $\delta_{C}^{em}$ and $\delta_{C}^{h}$, we can
employ the expressions which were derived in Ref. [11] to calculate
$\delta_{C}$.

###  3.1 Theoretical approach used [11] to calculate Coulomb corrections to
the non-flip amplitudes

Considering multiple photon exchange in the elastic pp scattering and
neglecting the higher order corrections ${\cal O}(\alpha^{3})$, the net long
range Coulomb (C) amplitude (see Fig.1) can be presented as [11]

$\displaystyle f_{C}(q_{T})$
$\displaystyle=\frac{i}{2\pi}\int{d^{2}b\,e^{i\vec{q}_{T}\vec{b}}\left[1-e^{i\chi_{C}^{\textit{nf}}(b)}\right]},$
(9)
$\displaystyle=\hat{f}_{C}(q_{T})+\frac{i}{2\pi}\int{d^{2}b\,e^{i\vec{q}_{T}\vec{b}}\left[\chi_{C}(b)\right]^{2}/2}$
(10)

where $q_{T}\\!\approx\\!\sqrt{-t}$ is transverse momentum and the eikonal
phase

$\chi_{C}^{\textit{nf}}(b)=\frac{1}{2\pi}\int{d^{2}q_{T}\hat{f}_{C}(q_{T})e^{-i\vec{q}_{T}\vec{b}}}$
(11)

is a Fourier transform of the Coulomb part of the amplitude calculated in Born
approximation [11]

$\hat{f}_{C}(q_{T})=\frac{-2\alpha}{q_{T}^{2}+\lambda^{2}}\,e^{-\widetilde{B}q_{T}^{2}/2}.$
(12)

Here, $\hat{f}_{C}$ is defined as sum of two non-flip helicity amplitudes
$\langle\\!+\>\\!\\!+\\!|\\!+\>\\!\\!+\\!\rangle$ and
$\langle\\!+\>\\!\\!-\\!|\\!+\>\\!\\!-\\!\rangle$ [15]. A small photon mass
$\lambda$ was included to Eq. (12) to keep the integrals finite.

The multi-photon exchange results in an acquired Coulomb phase
$\Phi_{C}(q_{T})$

$f_{C}(q_{T})=\hat{f}_{C}(q_{T})\,e^{i\Phi_{C}(q_{T})}.$ (13)

Assuming $\Phi_{C}(q_{T})\\!\ll\\!1$, one finds

$\displaystyle\\!\\!\Phi_{C}(q_{T})$
$\displaystyle=-i\left[f_{C}(q_{T})/\hat{f}_{C}(q_{T})-1\right]$ (14)
$\displaystyle=\frac{1}{4\pi}\int{d^{2}q_{1}d^{2}q_{2}\delta(\vec{q}_{T}\\!-\\!\vec{q}_{1}\\!-\\!\vec{q}_{2})\frac{\hat{f}_{C}(q_{1})\hat{f}_{C}(q_{2})}{\hat{f}_{C}(q_{T})}}\\!,$
(15)

Similarly, to calculate the Coulomb corrections to the hadronic amplitude,

$\hat{f}_{N}(q_{T})=\frac{(i+\rho)\sigma_{\text{tot}}}{4\pi}\,e^{-Bq_{T}^{2}/2},$
(16)

one can use the following relations:

$\displaystyle f_{\text{NC}}(q_{T})$
$\displaystyle=\frac{i}{2\pi}\int{d^{2}b\,e^{i\vec{q}_{T}\vec{b}}\gamma_{N}^{\textit{nf}}(b)e^{i\chi_{C}^{\textit{nf}}(b)}}$
(17)
$\displaystyle=\hat{f}_{N}(q_{T})+\frac{i}{2\pi}\int{d^{2}b\,e^{i\vec{q}_{T}\vec{b}}\gamma_{N}^{\textit{nf}}\chi_{C}^{\textit{nf}}},$
(18) $\displaystyle\gamma_{N}^{\textit{nf}}(b)$
$\displaystyle=\frac{-i}{2\pi}\int{d^{2}q_{T}e^{-i\vec{q}_{T}\vec{b}}\,\hat{f}_{N}(q_{T})},$
(19) $\displaystyle\Phi_{\text{NC}}(q_{T})$
$\displaystyle=-\frac{\alpha}{\pi}\int{\frac{d^{2}q_{1}}{q_{1}^{2}+\lambda^{2}}}$
$\displaystyle\times\exp{\left[-(B+\widetilde{B})q_{1}^{2}/2+B\vec{q}_{1}\vec{q_{T}}.\right]}$
(20)

Eqs. (15) and (20) were analytically integrated in Ref. [11]. Both,
$\Phi_{C}(q_{T})$ and $\Phi_{\text{NC}}(q_{T})$, contains the divergent term
$\ln{q^{2}/\lambda^{2}}$ which, however, cancels in final expression (4) for
the Coulomb phase difference

$\delta_{C}(t)=\Phi_{C}(t)-\Phi_{\text{NC}}(t).$ (21)

###  3.2 Coulomb corrections to the spin-flip amplitudes

To find the Coulomb corrected spin-flip amplitudes
$f_{C}^{\textit{sf}}(q_{T})$ and $f_{\text{NC}}^{\textit{sf}}(q_{T})$, one can
use the following eikonal phases [20]

$\chi_{C}^{\textit{sf}}(b)=\frac{1}{2\pi}\int{d^{2}qe^{-i\vec{q}\vec{b}}\times\frac{\kappa_{p}}{2m_{p}}\,(\vec{n}\vec{q})\,\hat{f}_{C}(q)/2}$
(22)

and

$\gamma_{N}^{\textit{sf}}(b)=\frac{-i}{2\pi}\int{d^{2}qe^{-i\vec{q}\vec{b}}\times\frac{r_{5}}{(i+\rho)m_{p}}\,(\vec{n}\vec{q})\,\hat{f}_{N}(q)/2}$
(23)

respectively. Here, $\vec{n}\propto\vec{p}_{\text{beam}}\times\vec{s}$ is a
unit vector orthogonal to the beam momentum $\vec{p}_{\text{beam}}$ and the
proton spin $\vec{s}$.

Considering the spin flip amplitudes for $\vec{q}_{T}\\!=\\!\vec{n}q_{T}$, one
can readily determine the spin-flip phase $\Phi_{C}^{\textit{sf}}(q_{T})$ by
including factor $2(\vec{q}_{T}\vec{q}_{1})/q_{T}^{2}$ or
$2(\vec{q}_{T}\vec{q}_{2})/q_{T}^{2}$ to integral (15). Since
$\vec{q}_{1}\vec{q}_{T}\\!+\\!\vec{q}_{2}\vec{q}_{T}\\!=\\!q_{T}^{2}$, we
immediately find

$\Phi_{C}^{\textit{sf}}(q_{T})\\!=\\!\Phi_{C}(q_{T}),$ (24)

which leads to

$\delta_{C}^{\textit{em}}(t)=\delta_{C}(t).$ (25)

To calculate $\Phi_{\text{NC}}^{\textit{sf}}(q_{T})$, factor
$(\vec{q}_{2}\vec{q}_{T})/q_{T}^{2}=1-(\vec{q}_{1}\vec{q}_{T})/q_{T}^{2}$
should be applied to Eq. (20), which gives

$\displaystyle\\!\\!\Phi_{\text{NC}}^{\textit{sf}}(q_{T})$
$\displaystyle=\Phi_{\text{NC}}(q_{T})-\frac{\alpha
B}{B+\widetilde{B}}\times\Delta_{\text{NC}}^{\textit{sf}}(\eta),$ (26)
$\displaystyle\Delta_{\text{NC}}^{\textit{sf}}(\eta)$
$\displaystyle=\int_{0}^{\infty}{\frac{du}{\eta}e^{-u^{2}/4\eta}}\int_{-\pi}^{\pi}{\frac{d\varphi}{2\pi}\cos{\varphi}\,e^{u\cos{\varphi}}},$
(27) $\displaystyle\eta$
$\displaystyle=\frac{B}{B+\widetilde{B}}\times\frac{Bq_{T}^{2}}{2}\approx
Bq_{T}^{2}/4.$ (28)

Expanding

$e^{u\cos{\varphi}}\cos{\varphi}\to\sum_{k=0}^{\infty}{\frac{u^{k}}{k!}\cos^{k+1}{\varphi}}$
(29)

and using following integrals [21]

$\displaystyle\int_{-\pi}^{\pi}{\cos^{2n+1}{x}\,dx}$ $\displaystyle=0,$ (30)
$\displaystyle\int_{-\pi}^{\pi}{\cos^{2n}{x}\,dx}$
$\displaystyle=\frac{\pi}{2^{n-2}}\,\frac{(2n-1)!}{(n-1)!n!},$ (31)
$\displaystyle\int_{0}^{\infty}{x^{2n+1}e^{-px^{2}}\,dx}$
$\displaystyle=\frac{n!}{2p^{n+1}},$ (32)

one arrives to

$\Delta_{\text{NC}}^{\textit{sf}}(\eta)=\sum_{k=0}^{\infty}{\frac{\eta^{k}}{(k+1)!}}=\frac{e^{\eta}-1}{\eta}.$
(33)

Thus,

$\delta_{C}^{h}(t)=\delta_{C}(t)-\frac{\alpha
B}{B+\widetilde{B}}\,\frac{e^{\eta}-1}{\eta}.$ (34)

In Ref. [7], it was pointed out that hadronic spin-flip amplitude should also
include the spin-flip photon exchange, i.e. one should replace

$\gamma_{N}^{\textit{sf}}(b)\to\gamma_{N}^{\textit{sf}}(b)+i\chi_{C}^{\textit{sf}}(b)\gamma_{N}^{\textit{nf}}(b)$
(35)

in Eq. (19). Using Eq. (26), one readily finds

$\displaystyle\frac{i}{2\pi}$
$\displaystyle\int{d^{2}be^{i\vec{q}_{T}\vec{b}}\,\chi_{C}^{\textit{sf}}(b)\gamma_{N}^{\textit{nf}}(b)}$
$\displaystyle=\frac{\kappa_{p}\,\left[\Phi_{N}(q_{T})\\!-\\!\Phi_{N}^{\textit{sf}}(q_{T})\right]}{2m_{p}}(\vec{n}\vec{q}_{T})\hat{f}_{N}(q_{T}).$
(36)

Thus, replacement (35) can be interpreted as an effective alteration of the
hadronic spin-flip parameter

$\displaystyle r_{5}\to$ $\displaystyle r_{5}+i(i+\rho)\Delta_{\gamma}\approx
r_{5}-\Delta_{\gamma},$ (37)
$\displaystyle\Delta_{\gamma}=\frac{\kappa_{p}}{2}\frac{\alpha
B}{(B+\widetilde{B})}\approx 0.003.$ (38)

For the modified $r_{5}$, Coulomb phase $\Phi_{\text{NC}}^{\textit{sf}}(t)$ is
the same as in Eq. (26).

##  4 Summary

The technique developed in Ref. [11] was adapted for calculation of the
Coulomb phase shifts in the spin-flip terms
$\widetilde{\phi}_{\textit{sf}}\phi_{\textit{nf}}^{*}$ (25) and
$\phi_{\textit{sf}}\widetilde{\phi}_{\textit{nf}}^{*}$ (34) of
$A_{\text{N}}(t)$.

Small difference, $\delta_{C}^{h}\\!-\\!\delta_{C}\sim-\alpha/2$, was found
for the CNI scattering, $|t|\lesssim 0.1\,\text{GeV}^{2}$ ($\eta<0.3$). Since
$|r_{5}|\\!\lesssim 0.02$, such a discrepancy can be neglected in Eq. (8).

Thus, we should agree with approximation for the Coulomb phases

$\delta_{C}^{\textit{em}}=\delta_{C}^{\textit{h}}=\delta_{C}=-\alpha\times\left[\ln{\frac{(B\\!+\\!\widetilde{B})|t|}{2}}+\gamma\right]$
(39)

suggested in Ref. [15].

Evaluating $\delta_{C}^{\textit{em}}$ and $\delta_{C}^{h}$, we did not
distinguish between non-flip $B$ and spin-flip $B_{\textit{sf}}$ hadronic
slopes as well as between $\widetilde{B}_{\textit{sf}}$ and

$\widetilde{B}_{\textit{sf}}=\left(r_{E}^{2}+r_{M}^{2}\right)/3,$ (40)

where $r_{M}=0.851\\!\pm\\!0.026\,\text{fm}$ [22] is rms magnetic radius of a
proton.

Using explicit expressions for $\Phi_{C}(t)$ and $\Phi_{N}(t)$ [11] one can
derive

$\displaystyle\\!\\!\delta_{C}^{\textit{em}}(t,\widetilde{B}_{\textit{sf}},B)$
$\displaystyle=\delta_{C}(t,\widetilde{B},B)+{\cal
O}(|\widetilde{B}_{\textit{sf}}\\!-\\!\widetilde{B}|t),$ (41)
$\displaystyle\delta_{C}^{h}(t,\widetilde{B},B_{\textit{sf}})$
$\displaystyle=\delta_{C}(t,\widetilde{B},B)-\frac{\alpha
B_{\textit{sf}}}{B_{\textit{sf}}+\widetilde{B}}$
$\displaystyle\quad+\alpha\ln{\frac{B+\widetilde{B}}{B_{\textit{sf}}+\widetilde{B}}}+{\cal
O}(|B_{\textit{sf}}\\!-\\!B|t).$ (42)

For a reasonable ratio $B_{\textit{sf}}/B$ between the spin-flip and non-flip
slopes [23], the related corrections to $\delta_{C}^{\textit{em}}$ and
$\delta_{C}^{h}$ can be neglected in the experimental data analysis.

Spin-flip photon exchange (35) results in an effective correction

$\Delta{\rm Re}\,r_{5}=-\Delta_{\gamma}\approx-0.003.$ (43)

to the value of the hadronic spin-flip parameter $r_{5}$ in Eq. (8). The
correction found is about triple of the experimental accuracy for $R_{5}$ in
the HJET measurements [16]. Since $\Delta_{\gamma}$ is independent of the
experimental data analysis, any already published result for $r_{5}$, in
particular [16, 19], can be easily adjusted

${\rm Re}\,r_{5}={\rm
Re}\,r_{5}^{\text{(publ.)}}+\frac{\alpha\kappa_{p}B}{2(B+\widetilde{B})}.$
(44)

This might be especially important for a study of $r_{5}(s)$ dependence on
energy, e.g. in the Regge fit [16, 7].

To illustrate a possible effect of corrections (44), we can reconsider the
Regge fit [16] of the HJET values of $r_{5}$, measured at
$\sqrt{s}\\!=\\!13.76~{}\text{and}~{}21.91\,\text{GeV}$:

$\sigma_{\text{tot}}(s)\times
r_{5}(s)=\sum_{{\cal{}R}=P,R^{\pm}}{f_{5}^{\cal{}R}\,{\cal{}R}(s)}.$ (45)

In Ref. [16], the spin-flip Reggeon, $R^{+}(s)$ and $R^{-}(s)$, and Pomeron
(in a Froissaron parametrization), $P(s)$, functions were approximated by the
non-flip ones [3]. After applying corrections (44), the fit
$\chi^{2}\\!=\\!2.2$ (NDF=1) was improved to $\chi^{2}\\!=\\!0.9$ and the
central values of couplings $f_{5}^{\cal{}R}$ were shifted to

$\displaystyle\phantom{-}f_{5}^{P^{\phantom{+}}}\\!\\!\\!=\\!0.045\\!\pm\\!0.002_{\text{stat}}\\!\pm\\!0.003_{\text{syst}}$
$\displaystyle\,\to\,0.053,$ (46) $\displaystyle-
f_{5}^{R^{+}}\\!\\!\\!=\\!0.032\\!\pm\\!0.007_{\text{stat}}\\!\pm\\!0.014_{\text{syst}}$
$\displaystyle\,\to\,0.069,$ (47)
$\displaystyle\phantom{-}f_{5}^{R^{-}}\\!\\!\\!=\\!0.622\\!\pm\\!0.023_{\text{stat}}\\!\pm\\!0.024_{\text{syst}}$
$\displaystyle\,\to\,0.654.$ (48)

A thorough re-fit of the $r_{5}$ measurements, including the STAR
$\sqrt{s}\\!=\\!200\,\text{GeV}$ result [19], will be done elsewhere.

It is interesting to note that an absorptive correction,

$a_{\textit{sf}}=\alpha B/(B+\widetilde{B}_{\textit{sf}}),$ (49)

to the electromagnetic spin-flip form factor, ${\cal
F}_{pp}^{\textit{sf}}(t)\to{\cal
F}_{pp}^{\textit{sf}}(t)\times\left(1+a_{\textit{sf}}\,t/t_{c}\right)$,
results [6] in the same effective alteration of $r_{5}$ as shown in Eq. (43).

In Ref. [7], to calculate absorptive corrections, graphs in Fig. 1 were re-
grouped. In this approach, a spin-flip photon contribution [Eq. (35] to
hadronic spin-flip amplitude $f_{N}^{\textit{sf}}$ was not considered.
However, the term $\chi_{C}^{\textit{sf}}(b)\gamma_{N}^{\textit{nf}}(b)$
appeared as an absorptive correction to electromagnetic spin-flip amplitude
$f_{C}^{\textit{sf}}(q_{T})$. The corresponding equation was integrated
numerically in [7]. Here, using (36), we found absorptive correction parameter
$a_{\textit{sf}}$ to be identical to that of given in Eq. (49).

## Acknowledgements

The author is thankful to Nigel Buttimore and Boris Kopeliovich for
informative discussions and acknowledges support from the Office of Nuclear
Physics in the Office of Science of the US Department of Energy.

## References

  * [1] A. I. Akhiezer, I. Y. Pomeranchuk, About elastic scattering of fast charged particles by nuclei, Zh. Eksp. Teor. Fiz. 16 (1946) 396, [J. Phys. (USSR) 9, 471 (1945)].
  * [2] B. Kopeliovich, L. Lapidus, On the necessity of polarization experiments in colliding p p and anti-p p beams, Yad. Fiz. 19 (1974) 218–223, [Sov. J. Nucl. Phys. 19, 114 (1974)].
  * [3] D. Fagundes, M. Menon, P. Silva, Leading components in forward elastic hadron scattering: Derivative dispersion relations and asymptotic uniqueness, Int. J. Mod. Phys. A 32 (32) (2017) 1750184.  arXiv:1705.01504, doi:10.1142/S0217751X17501846.
  * [4] V. Bartenev, et al., Measurement of the Slope of the Diffraction Peak for Elastic pp Scattering from 8-GeV to 400-GeV., Phys. Rev. Lett. 31 (1973) 1088–1091.  doi:10.1103/PhysRevLett.31.1088.
  * [5] V. Bartenev, et al., Real Part of the Proton Proton Forward Scattering Amplitude from 50 GeV to 400 GeV., Phys. Rev. Lett. 31 (1973) 1367–1370.  doi:10.1103/PhysRevLett.31.1367.
  * [6] A. A. Poblaguev, Corrections to the Elastic Proton-Proton Analyzing Power Parametrization at High Energies, Phys. Rev. D 100 (11) (2019) 116017.  arXiv:1910.02563, doi:10.1103/PhysRevD.100.116017.
  * [7] B. Z. Kopeliovich, M. Krelina, I. K. Potashnikova, Probing the Pomeron spin structure with Coulomb-nuclear interference, Phys. Lett. B 816 (2021) 136262\.  arXiv:2102.01595, doi:10.1016/j.physletb.2021.136262.
  * [8] H. A. Bethe, Scattering and polarization of protons by nuclei, Annals Phys. 3 (1958) 190–240.  doi:10.1016/0003-4916(58)90017-4.
  * [9] G. B. West, D. R. Yennie, Coulomb interference in high-energy scattering, Phys. Rev. 172 (1968) 1413–1422.  doi:10.1103/PhysRev.172.1413.
  * [10] R. Cahn, Coulombic - Hadronic Interference in an Eikonal Model, Z. Phys. C 15 (1982) 253.  doi:10.1007/BF01475009.
  * [11] B. Z. Kopeliovich, A. V. Tarasov, The Coulomb phase revisited, Phys. Lett. B 497 (2001) 44–48.  arXiv:hep-ph/0010062, doi:10.1016/S0370-2693(00)01316-2.
  * [12] L. H. Chan, K. W. Chen, J. R. Dunning, N. F. Ramsey, J. K. Walker, R. Wilson, Nucleon Form Factors and Their Interpretation, Phys. Rev. 141 (1966) 1298–1307.  doi:10.1103/PhysRev.141.1298.
  * [13] P. Zyla, et al., Review of Particle Physics, PTEP 2020 (8) (2020) 083C01.  doi:10.1093/ptep/ptaa104.
  * [14] N. H. Buttimore, E. Gotsman, E. Leader, Spin-dependent phenomena induced by electromagnetic-hadronic interference at high energies, Phys. Rev. D 18 (1978) 694–716, [Erratum: 35, 407 (1987)].  doi:10.1103/PhysRevD.18.694.
  * [15] N. H. Buttimore, B. Kopeliovich, E. Leader, J. Soffer, T. Trueman, The spin dependence of high-energy proton scattering, Phys. Rev. D 59 (1999) 114010.  arXiv:hep-ph/9901339, doi:10.1103/PhysRevD.59.114010.
  * [16] A. A. Poblaguev, et al., Precision Small Scattering Angle Measurements of Elastic Proton-Proton Single and Double Spin Analyzing Powers at the RHIC Hydrogen Jet Polarimeter, Phys. Rev. Lett. 123 (16) (2019) 162001.  arXiv:1909.11135, doi:10.1103/PhysRevLett.123.162001.
  * [17] H. Okada, et al., Measurement of the analyzing power in pp elastic scattering in the peak CNI region at RHIC, Phys. Lett. B 638 (2006) 450–454.  arXiv:nucl-ex/0502022, doi:10.1016/j.physletb.2006.06.008.
  * [18] I. G. Alekseev, et al., Measurements of single and double spin asymmetry in pp elastic scattering in the CNI region with a polarized atomic hydrogen gas jet target, Phys. Rev. D 79 (2009) 094014.  doi:10.1103/PhysRevD.79.094014.
  * [19] L. Adamczyk, et al., Single Spin Asymmetry $A_{N}$ in Polarized Proton-Proton Elastic Scattering at $\sqrt{s}=200$ GeV, Phys. Lett. B 719 (2013) 62–69.  arXiv:1206.1928, doi:10.1016/j.physletb.2013.01.014.
  * [20] L. I. Lapidus, Polarization Phenomena in Hadronic Collisions at Small Momentum Transfer, Fiz. Elem. Chast. Atom. Yadra 9 (1978) 84–117.
  * [21] I. S. Gradshteyn, I. M. Ryzhik, D. Zwillinger, V. Moll, Table of integrals, series, and products; 8th ed., Academic Press, Amsterdam, 2014.  doi:10.1016/C2010-0-64839-5.
  * [22] G. Lee, J. R. Arrington, R. J. Hill, Extraction of the proton radius from electron-proton scattering data, Phys. Rev. D 92 (1) (2015) 013013.  arXiv:1505.01489, doi:10.1103/PhysRevD.92.013013.
  * [23] J. R. Cudell, E. Predazzi, O. V. Selyugin, High-energy hadron spin flip amplitude at small momentum transfer and new A(N) data from RHIC, Eur. Phys. J. A 21 (2004) 479–486.  arXiv:hep-ph/0401040, doi:10.1140/epja/i2004-10012-2.


* (62) W. Chen, Reduction of Feynman Integrals in the Parametric Representation, JHEP 02 (2020) 115, [arXiv:1902.10387].
  * (63) W. Chen, Reduction of Feynman Integrals in the Parametric Representation II: Reduction of Tensor Integrals, Eur. Phys. J. C 81 (2021), no. 3 244, [arXiv:1912.08606].
  * (64) W. Chen, Reduction of Feynman integrals in the parametric representation III: integrals with cuts, Eur. Phys. J. C 80 (2020), no. 12 1173, [arXiv:2007.00507].
  * (65) T. D. Lee and M. Nauenberg, Degenerate systems and mass singularities, Phys. Rev. 133 (Mar, 1964) B1549–B1562.
  * (66) T. Kinoshita, Mass singularities of Feynman amplitudes, J. Math. Phys. 3 (1962) 650–677.
  * (67) R. N. Lee, Presenting LiteRed: a tool for the Loop InTEgrals REDuction, arXiv:1212.2685.
  * (68) R. N. Lee, LiteRed 1.4: a powerful tool for reduction of multiloop integrals, J. Phys. Conf. Ser. 523 (2014) 012059, [arXiv:1310.1145].
  * (69) A. V. Smirnov and F. S. Chuharev, FIRE6: Feynman Integral REduction with Modular Arithmetic, Comput. Phys. Commun. 247 (2020) 106877, [arXiv:1901.07808].
  * (70) C. Meyer, Algorithmic transformation of multi-loop master integrals to a canonical basis with CANONICA, Comput. Phys. Commun. 222 (2018) 295–312, [arXiv:1705.06252].
  * (71) R. N. Lee, Reducing differential equations for multiloop master integrals, JHEP 04 (2015) 108, [arXiv:1411.0911].
  * (72) R. N. Lee, Libra: A package for transformation of differential systems for multiloop integrals, Comput. Phys. Commun. 267 (2021) 108058, [arXiv:2012.00279].
  * (73) J. M. Campbell and E. W. N. Glover, Double unresolved approximations to multiparton scattering amplitudes, Nucl. Phys. B 527 (1998) 264–288, [hep-ph/9710255].
  * (74) S. Catani and M. Grazzini, Collinear factorization and splitting functions for next-to-next-to-leading order QCD calculations, Phys. Lett. B 446 (1999) 143–152, [hep-ph/9810389].
  * (75) M. Ritzmann and W. J. Waalewijn, Fragmentation in Jets at NNLO, Phys. Rev. D 90 (2014), no. 5 054029, [arXiv:1407.3272].
  * (76) A. Gehrmann-De Ridder and E. W. N. Glover, A Complete O (alpha alpha-s) calculation of the photon + 1 jet rate in e+ e- annihilation, Nucl. Phys. B 517 (1998) 269–323, [hep-ph/9707224].
  * (77) F. V. Tkachov, A Theorem on Analytical Calculability of Four Loop Renormalization Group Functions, Phys. Lett. B 100 (1981) 65–68.
  * (78) K. G. Chetyrkin and F. V. Tkachov, Integration by Parts: The Algorithm to Calculate beta Functions in 4 Loops, Nucl. Phys. B 192 (1981) 159–204.
  * (79) R. N. Lee, Modern techniques of multiloop calculations, in 49th Rencontres de Moriond on QCD and High Energy Interactions, pp. 297–300, 2014\.  arXiv:1405.5616.
  * (80) P. Maierhöfer, J. Usovitsch, and P. Uwer, Kira—A Feynman integral reduction program, Comput. Phys. Commun. 230 (2018) 99–112, [arXiv:1705.05610].
  * (81) J. Klappert and F. Lange, Reconstructing rational functions with FireFly, Comput. Phys. Commun. 247 (2020) 106951, [arXiv:1904.00009].
  * (82) J. Klappert, S. Y. Klein, and F. Lange, Interpolation of dense and sparse rational functions and other improvements in FireFly, Comput. Phys. Commun. 264 (2021) 107968, [arXiv:2004.01463].
  * (83) J. Klappert, F. Lange, P. Maierhöfer, and J. Usovitsch, Integral reduction with Kira 2.0 and finite field methods, Comput. Phys. Commun. 266 (2021) 108024, [arXiv:2008.06494].
  * (84) A. V. Kotikov, Differential equations method: New technique for massive Feynman diagrams calculation, Phys. Lett. B 254 (1991) 158–164.
  * (85) E. Remiddi, Differential equations for Feynman graph amplitudes, Nuovo Cim. A 110 (1997) 1435–1452, [hep-th/9711188].
  * (86) J. M. Henn, Multiloop integrals in dimensional regularization made simple, Phys. Rev. Lett. 110 (2013) 251601, [arXiv:1304.1806].
  * (87) A. B. Goncharov, Multiple polylogarithms and mixed tate motives, 2001.
  * (88) A. B. Goncharov, Multiple polylogarithms, cyclotomy and modular complexes, Math. Res. Lett. 5 (1998) 497–516, [arXiv:1105.2076].
  * (89) J. M. Borwein, D. M. Bradley, D. J. Broadhurst, and P. Lisonek, Special values of multiple polylogarithms, Trans. Am. Math. Soc. 353 (2001) 907–941, [math/9910045].
  * (90) L. Naterop, A. Signer, and Y. Ulrich, handyG —Rapid numerical evaluation of generalised polylogarithms in Fortran, Comput. Phys. Commun. 253 (2020) 107165, [arXiv:1909.01656].
  * (91) R. E. Cutkosky, Singularities and discontinuities of Feynman amplitudes, J. Math. Phys. 1 (1960) 429–433.
  * (92) C. Anastasiou and K. Melnikov, Higgs boson production at hadron colliders in NNLO QCD, Nucl. Phys. B 646 (2002) 220–256, [hep-ph/0207004].
  * (93) D. Maitre, HPL, a mathematica implementation of the harmonic polylogarithms, Comput. Phys. Commun. 174 (2006) 222–240, [hep-ph/0507152].
  * (94) A. Gehrmann-De Ridder, T. Gehrmann, and G. Heinrich, Four particle phase space integrals in massless QCD, Nucl. Phys. B 682 (2004) 265–288, [hep-ph/0311276].
  * (95) V. Magerya and A. Pikelner, Cutting massless four-loop propagators, JHEP 12 (2019) 026, [arXiv:1910.07522].
  * (96) G. P. Korchemsky and G. F. Sterman, Power corrections to event shapes and factorization, Nucl. Phys. B 555 (1999) 335–351, [hep-ph/9902341].
  * (97) Y. L. Dokshitzer, G. Marchesini, and B. R. Webber, Nonperturbative effects in the energy energy correlation, JHEP 07 (1999) 012, [hep-ph/9905339].
  * (98) S. T. Schindler, I. W. Stewart, and Z. Sun, Renormalons in the energy-energy correlator, arXiv:2305.19311.
  * (99) T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen, and P. Z. Skands, An introduction to PYTHIA 8.2, Comput. Phys. Commun. 191 (2015) 159–177, [arXiv:1410.3012].
  * (100) Z. Tulipánt, A. Kardos, and G. Somogyi, Energy–energy correlation in electron–positron annihilation at NNLL + NNLO accuracy, Eur. Phys. J. C 77 (2017), no. 11 749, [arXiv:1708.04093].
  * (101) A. Kardos, S. Kluth, G. Somogyi, Z. Tulipánt, and A. Verbytskyi, Precise determination of $\alpha_{S}(M_{Z})$ from a global fit of energy–energy correlation to NNLO+NNLL predictions, Eur. Phys. J. C 78 (2018), no. 6 498, [arXiv:1804.09146].
  * (102) H. Cao, X. Liu, and H. X. Zhu, Toward precision measurements of nucleon energy correlators in lepton-nucleon collisions, Phys. Rev. D 107 (2023), no. 11 114008, [arXiv:2303.01530].
  * (103) O. V. Tarasov, A. A. Vladimirov, and A. Y. Zharkov, The Gell-Mann-Low Function of QCD in the Three Loop Approximation, Phys. Lett. B 93 (1980) 429–432.
  * (104) S. A. Larin and J. A. M. Vermaseren, The Three loop QCD Beta function and anomalous dimensions, Phys. Lett. B 303 (1993) 334–336, [hep-ph/9302208].
  * (105) T. van Ritbergen, J. A. M. Vermaseren, and S. A. Larin, The Four loop beta function in quantum chromodynamics, Phys. Lett. B 400 (1997) 379–384, [hep-ph/9701390].
  * (106) M. Czakon, The Four-loop QCD beta-function and anomalous dimensions, Nucl. Phys. B 710 (2005) 485–498, [hep-ph/0411261].
  * (107) T. Gehrmann, E. W. N. Glover, T. Huber, N. Ikizlerli, and C. Studerus, Calculation of the quark and gluon form factors to three loops in QCD, JHEP 06 (2010) 094, [arXiv:1004.3653].