# Fermionic Sign Problem Minimization by Constant Path Integral Contour Shifts

Christoph Gäntgen  Institute for Advanced Simulation, Forschungszentrum
Jülich, 54245 Jülich, Germany  Center for Advanced Simulation and Analytics
(CASA), Forschungszentrum Jülich, 52425 Jülich, Germany  Helmholtz-Institut
für Strahlen- und Kernphysik, Rheinische Friedrich-Wilhelms-Universität Bonn,
53115 Bonn, Germany  Evan Berkowitz  Institute for Advanced Simulation,
Forschungszentrum Jülich, 54245 Jülich, Germany  Center for Advanced
Simulation and Analytics (CASA), Forschungszentrum Jülich, 52425 Jülich,
Germany  JARA & Jülich Supercomputing Center, Forschungszentrum Jülich, 54245
Jülich, Germany  Thomas Luu  Institute for Advanced Simulation,
Forschungszentrum Jülich, 54245 Jülich, Germany  Helmholtz-Institut für
Strahlen- und Kernphysik, Rheinische Friedrich-Wilhelms-Universität Bonn,
53115 Bonn, Germany  Institut für Kernphysik, Forschungszentrum Jülich, 54245
Jülich, Germany  Johann Ostmeyer  Department of Mathematical Sciences,
University of Liverpool, Liverpool, L69 7ZL, United Kingdom  Marcel Rodekamp
Institute for Advanced Simulation, Forschungszentrum Jülich, 54245 Jülich,
Germany  Center for Advanced Simulation and Analytics (CASA),
Forschungszentrum Jülich, 52425 Jülich, Germany  Helmholtz-Institut für
Strahlen- und Kernphysik, Rheinische Friedrich-Wilhelms-Universität Bonn,
53115 Bonn, Germany  JARA & Jülich Supercomputing Center, Forschungszentrum
Jülich, 54245 Jülich, Germany

###### Abstract

The path integral formulation of quantum mechanical problems including
fermions is often affected by a severe numerical sign problem. We show how
such a sign problem can be alleviated by a judiciously chosen constant
imaginary offset to the path integral. Such integration contour deformations
introduce no additional computational cost to the Hybrid Monte Carlo
algorithm, while its effective sample size is greatly increased. This makes
otherwise unviable simulations efficient for a wide range of parameters.
Applying our method to the Hubbard model, we find that the sign problem is
significantly reduced. Furthermore, we prove that it vanishes completely for
large chemical potentials, a regime where the sign problem is expected to be
particularly severe without imaginary offsets. In addition to a numerical
analysis of such optimized contour shifts, we analytically compute the shifts
corresponding to the leading and next-to-leading order corrections to the
action. We find that such simple approximations, free of significant
computational cost, suffice in many cases.

##  I Introduction

The _numerical sign problem_ is a major hindrance for the application of
stochastic methods to certain physical systems, such as QCD at finite baryon
density or electronically doped systems in strongly correlated condensed
matter. The problem refers to the extreme cost of numerically approximating
integrals arising with a highly oscillatory integrand, such as path integrals
with complex-valued actions. Because partition functions are exponential in
the action, the numerical costs typically scale exponentially in the spacetime
volume [1], pushing many physically interesting systems beyond the reach of
numerical investigation.

Methods that reduce the sign problem allow us to use our limited resources
more efficiently and thus extend the range of systems we can investigate. In
cases when the offending term in the Hamiltonian that induces the complex
phase is small, one can rely on simple reweighting. For small systems or
ground-state properties one can forego stochastic simulations and instead use
direct methods, such as tensor networks [2]. Complex Langevin is another
popular method to fight the sign problem, but a lot of technology is required
to guarantee it converges to the right distribution [3, 4]. Each of these
methods have their own limitations, by no means fully solving the sign
problem.

Here we focus on _contour deformation_ to alleviate the sign problem. One
transforms the integration domain of the path integral to a more favorable
manifold in the high-dimensional complex space where the sign oscillations are
reduced [5, 6, 7, 8, 9, 10, 11, 4, 12, 13]. Such deformations are formally
allowed as long as one does not cross any singularities of the integrand and
one preserves the homology class of the integral. There exist manifolds, so-
called _Lefschetz thimbles_ , where the complex phase remains fixed. In
theories where one thimble dominates, the sign problem is solved since the
constant complex phase on the thimble can be factored outside of the path
integral. Even when multiple thimbles contribute, each with a different but
constant phase, the sign problem is not eliminated, but is expected to be
improved. While the locations of these thimbles are not known _a priori_ they
can be found by integrating holomorphic flow equations. Unfortunately the
numerical determination of the complete set of contributing thimbles is quite
costly: mapping out their full constellation is just as difficult as the
original sign problem.

As the goal in our work is to _alleviate_ the sign problem (as opposed to
eliminating it), a natural question arises: given finite computational
resources, which contour deformations are most efficient for the problem at
hand at alleviating the sign problem _sufficiently_ , meaning that observables
can be extracted in a statistically meaningful and reliable manner. In
previous studies of the Hubbard model [14, 13] we trained neural networks
(NNs) on flowed configurations to parametrize an integration manifold called a
_learnifold_ [15, 11]. This deformation worked very well at alleviating the
sign problem for various doped Hubbard systems. In particular, we provided
physical results of a doped Hubbard model for carbon nano-systems up to 18 ion
sites [13]. However this method still comes at the cost of generating flowed
training data and training a neural network. Here we study the simplest
imaginable deformation: shifting the integration manifold by a global
imaginary constant offset. We find that optimizing the offset substantially
reduces the computational demands and often yields a contour deformation of
equal potency.

A constant offset induces no Jacobian; keeping the method simple. Further, a
constant offset does not require modification of the Monte Carlo algorithm,
nor does it require generation of training data and training of NNs. In ref.
[15], for example, it was shown for the Thirring model that a calculation on
the _tangent plane_ , a constant offset that intersects the classical saddle
point of the main Lefschetz thimble, is sufficient at alleviating the sign
problem. We have also used this deformation as a comparison to our neural
networks in previous publications [14, 13]. For the Hubbard model, however, we
find that for certain values of the chemical potential, the tangent plane does
not meaningfully alleviate the sign problem. We will show how to incorporate
quantum corrections to the saddle point, thereby obtaining a better constant
shift that corresponds to the effective action obtained by inclusion of
1-particle irreducible terms. Even then there are cases where we resort to
numerical optimization of the imaginary offset.

We study the fermionic Hubbard model. As opposed to the systems investigated
in our earlier work [13], here we also consider systems that are non-
bipartite, such as the fullerenes $C_{20}$ and $C_{60}$. In the following
section we provide the formal aspects of our method, providing derivations for
the location of the classical and quantum-corrected saddle points that we use
to determine our constant offsets. In section III we continue with the
description of our method for numerically determining the optimized plane. We
then demonstrate the efficacy of our methods by providing numerical results of
various Hubbard systems in section IV. We recapitulate in section V. To keep
the presentation reasonable, we place formal (and tedious) derivations in the
appendices.

##  II Formalism

###  II.1 The Hubbard Model

The Hubbard model describes the interacting behavior of particles on a
lattice. In our case these are electrons on a lattice of ions. It consists of
a tight binding term and an onsite interaction representing electron-electron
repulsion. It takes into account external influences on the overall particle
number, like doping or an applied voltage, via a chemical potential $\mu$[16,
17, 18, 19, 20, 21, 22]. We formulate our theory in the _particle-hole basis_
[23]

$\displaystyle H$ $\displaystyle=-{\kappa}\sum_{\langle
x,y\rangle}\left(a^{\dagger}_{x}a_{y}-b^{\dagger}_{x}b_{y}\right)+\frac{U}{2}\sum_{x}q_{x}^{2}-{\mu}\sum_{x}q_{x}$
$\displaystyle q_{x}$
$\displaystyle=a{}^{\dagger}_{x}a_{x}-b{}^{\dagger}_{x}b_{x}$ (1)

where $\kappa$ is the hopping parameter of neighboring lattice sites, $U$
provides the strength of interaction of two electrons sharing one lattice
site, and $q$ is the local charge operator relative to half filling.
Alternatively the sum over neighboring sites can be represented with the
hopping matrix $K$ which is $\kappa$ times the adjacency matrix of the
lattice. This option also allows for individual hopping parameters. The sum in
$x$ is over all $N_{x}$ sites. The $a$ $(a^{\dagger})$ operator implements
particle destruction (creation) and $a^{\dagger}a$ counts particles. Similarly
the $b$ $(b^{\dagger})$ operators destroy (create) holes and $b^{\dagger}b$
counts them.

###  II.2 Finite lattices considered in this work

As we explain below, stochastic simulations of the Hubbard model suffer from
the sign problem when the geometry of the system is non-bipartite and/or a
non-zero chemical potential is present. A lattice is _bipartite_ when its
sites can be divided into two groups, such that each site has only neighbors
of the other group. Another way to think of it is that each closed path must
traverse an even number of links. In this paper we will investigate both
cases, with the 8- and 18-site honeycomb lattices as bipartite examples and
the $C_{20}$ and $C_{60}$ fullerenes as non-bipartite examples. The 8- and
18-site honeycomb lattices consist of $2\times 2$ and $3\times 3$ unit cells
respectively and in this work are assumed to have periodic boundary
conditions. $C_{20}$ is a dodecahedron with 12 equal pentagons. $C_{60}$ is a
truncated icosahedron with 12 pentagons and 20 hexagons. The four lattice
structures are visualised in fig. 1. All of the lattices we consider are _site
transitive_ , meaning that symmetries of the lattice can map any site to any
other site, an analog to translation invariance.

(a) 8-site

(b) 18-site

(c) $C_{20}$

(d) $C_{60}$

Figure 1: Spatial lattices considered in this paper.

###  II.3 The path-integral formulation of the Hubbard model

The expectation value of any quantum mechanical operator $\hat{O}$ can be
calculated within the path integral formalism,

$\expectationvalue{\hat{O}}=\frac{1}{\mathcal{Z}}\int\mathcal{D}\phi\,\hat{O}\left[\phi\right]e^{-S\left[\phi\right]}$
(2)

where the partition function is

$\mathcal{Z}=\int\mathcal{D}\phi\,e^{-S\left[\phi\right]}\ $ (3)

and
$\mathcal{D}\phi=\lim_{N_{t}\rightarrow\infty}\prod_{x}^{N_{x}}\prod_{t}^{N_{t}}d\phi_{x,t}$
and $S[\phi]$ is the action that defines the system. It is common practice to
estimate expectation values (2) with importance sampling

$\int\mathcal{D}\phi\,\hat{O}\left[\phi\right]\mathbb{P}\left[\phi\right]=\lim_{N\rightarrow\infty}\frac{1}{N}\sum_{i}^{N}\hat{O}\left[\phi_{i}\right],\quad\phi_{i}\sim\mathbb{P}\left[\phi_{i}\right],$
(4)

drawing configurations $\phi$ from the probability distribution $\mathbb{P}$
by Markov Chain Monte Carlo (MCMC) methods such as Hybrid (or Hamiltonian)
Monte Carlo (HMC) [24]. For purely real actions it is straight forward to
choose $\mathbb{P}\left[\phi\right]=e^{-S\left[\phi\right]}/\mathcal{Z}$. For
generally complex actions $S[\phi]=S_{R}[\phi]+iS_{I}[\phi]$, however, one
typically separates the complex phase and absorbs it into the definition of
the observable,

$\expectationvalue{\hat{O}}=\frac{\expectationvalue{\hat{O}\exp(-iS_{I})}_{R}}{\expectationvalue{\exp(-iS_{I})}_{R}}$
(5)

where the subscript _R_ indicates sampling according to the probability
defined by the real part of the action,
$\mathbb{P}[\phi]=e^{-S_{R}[\phi]}/\mathcal{Z}_{R}$. This process is called
_reweighting_ and it exactly produces correct expectation values (2) in the
limit of infinite statistics.

However, for finite statistics, the oscillating phase in the denominator of
the reweighting (5), known as the _average phase_ ,

$\expectationvalue{\exp(-iS_{I})}_{R}=\dfrac{\int\mathcal{D}\phi\,e^{-iS_{I}\left[\phi\right]}e^{-S_{R}\left[\phi\right]}}{\int\mathcal{D}\phi\,e^{-S_{R}\left[\phi\right]}}=\frac{\mathcal{Z}}{\mathcal{Z}_{R}}$
(6)

can be hard to numerically estimate. Stronger oscillation and its attendant
cancellations become more severe for larger system sizes and some parameters.
We call the average phase’s absolute value

$\Sigma=\absolutevalue{\expectationvalue{\exp(-iS_{I})}_{R}}=\absolutevalue{\frac{\mathcal{Z}}{\mathcal{Z}_{R}}}$
(7)

the _statistical power_ and we use it to quantify the sign problem. When the
statistical power is 1, the path integral is sign-problem free; a value of 0
indicates the worst possible sign problem. While in a problem-free case the
stochastic uncertainties of expectation values from an ensemble of
$N_{\text{cfg}}$ generated configurations scale with $N^{-1/2}_{\text{cfg}}$,
with a sign problem the statistical power effectively reduces the contribution
of each configuration: the error scales with the square root of the effective
number of samples [4],

$N_{\text{eff}}=\Sigma^{2}\times N_{\text{cfg}}\ .$ (8)

Trotterizing the thermal partition function $\tr{e^{-\beta H}}$, linearizing
the interaction with a Hubbard-Stratonovich transformation with an integral
over auxiliary fields $\phi$, and inserting resolutions of the identity in
terms of Grassmann coherent states, yields an action

$S\left[\phi,\tilde{K},\tilde{\mu}\right]=\frac{1}{2\tilde{U}}\sum_{t,x}\phi^{2}_{x,t}-\log\det\left(M\left[+\phi,+\tilde{K},+\tilde{\mu}\right]\right)-\log\det\left(M\left[-\phi,-\tilde{K},-\tilde{\mu}\right]\right)$
(9)

where the dimensionless parameters $\tilde{U}=U\times\delta$,
$\tilde{\kappa}=\kappa\delta$ etc. are rescaled by the temporal lattice
spacing $\delta=\beta/N_{t}$. The fermion matrices $M$ encode the particle and
hole fermion loops exactly. There are a variety discretizations of the fermion
matrix that become exact and equal in the continuum limit [9, 25, 23]. We use
the ‘exponential’ discretization in the language of ref. [9]

$M_{x^{\prime}t^{\prime},xt}[\pm\phi,\pm K,\pm\mu]\equiv
M_{x^{\prime}t^{\prime},xt}[\pm;\phi]=\delta_{x^{\prime},x}\delta_{t^{\prime},t}-[\exp(\pm\tilde{K})]_{x^{\prime},x}e^{\pm(-i\phi_{x,t}+\tilde{\mu})}B_{t^{\prime}}\delta_{t^{\prime},t+1}$
(10)

where $B_{N_{t}-1}$ carries an extra $-1$ encoding the fermionic temporal
antiperiodic boundary conditions. We do not consider other discretizations in
this work. We note, however, that when $\mu\neq 0$, this discretization does
not suffer from ergodicity issues described in refs. [25, 9].

We now consider deforming our original integral by complexifying the auxiliary
field $\phi$ and deforming the integration manifold. Cauchy’s theorem
guarantees that this deformation leaves all holomorphic observables the same,
as long as the deformation does not cross any singularities and the
deformation preserves the homology class. The statistical power depends on the
imaginary part of the action weighted by its real part and thus is not
holomorphic, so it is manifold-dependent. Some manifolds may tame the
oscillations, especially when they resemble Lefschetz thimbles, high-
dimensional analogs of contours of steepest descent [26, 11]. We have
previously trained neural networks to learn the results of the holomorphic
flow in a computationally tractable way [14, 13]. An even simpler deformation,
that of a constant imaginary shift in all components of $\phi$, can lead to
significant alleviation of the sign problem while incurring no additional
costs to the HMC algorithm. In particular, ref. [11] showed that a constant
shift that intersected the saddle point of the main thimble, producing the so
called ‘tangent plane’, sufficiently reduced the sign problem in simulations
of the Thirring model. We now consider the same constant shift to the tangent
plane of the Hubbard model.

###  II.4 The tangent plane of the Hubbard model

The holomorphic flow of a configuration $\phi$ is its image under evolution in
a fictitous time $t$ by

$\displaystyle\frac{d\phi}{dt}=(\partial_{\phi}S)^{*}.$ (11)

The saddle point that fixes the tangent plane is found by flowing the $\phi=0$
translationally-invariant configuration to its fixed point. Because it is a
fixed point the time derivative vanishes and the saddle point $\phi_{c}$
satisfies

$\displaystyle\left.\partial_{\phi}S[\phi]\right.\vline_{\phi=\phi_{c}}=0\ .$
(12)

In the graphene case this saddle point has the greatest weight [27] on the
semimetal side of the quantum critical point at $U\lesssim$3.8$$ [28, 29, 30,
31, 32, 33, 34]. The saddle point $\phi_{c}$ has zero real part and, because
the lattices we consider are site-transitive, constant imaginary part which is
non-zero when $\mu\neq 0$ on bipartite lattices and generically on non-
bipartite lattices.

Leveraging the simplicity of $\phi_{c}=i\phi_{0}$ independent of space and
time we can calculate the action

$\displaystyle
S[\phi_{c}=i\phi_{0}]=\frac{1}{2\tilde{U}}N_{x}N_{t}\phi_{0}^{2}-\log\det\left(\mathbbm{1}+e^{N_{t}\phi_{0}+\beta\mu}e^{+\beta
K}\right)-\log\det\left(\mathbbm{1}+e^{-N_{t}\phi_{0}-\beta\mu}e^{-\beta
K}\right)$ (13)

where we used the Schur complement to simplify the fermion determinants and
used the spatial independence of $\phi_{0}$ to commute the auxiliary field
terms past the hopping terms. Using $\log\det=\tr\log$, transforming to the
basis where the hopping matrix is diagonalized, and requiring $\phi_{c}$ to be
a fixed point (12) leads to

$\displaystyle\phi_{0}/\delta=-\frac{U}{N_{x}}\sum_{k}\tanh\left(\frac{\beta}{2}\left[\epsilon_{k}+\mu+\phi_{0}/\delta\right]\right)$
(14)

where the sum is over the $N_{x}$ modes of the hopping matrix $K$ with
noninteracting energy eigenvalue $\epsilon_{k}$; we provide a detailed
derivation in Appendix B. Writing $\delta=\beta/N_{t}$ shows that this
transcendental equation contains only temporal continuum quantities, except
for the combination $\phi_{0}N_{t}$, which will stay fixed as we go towards
the time continuum limit $N_{t}\rightarrow\infty$. We see that
$\phi_{0}/\tilde{U}$ is bounded between $-1$ and $+1$ and can cheaply
determine the imaginary offset of the tangent plane $\phi_{0}$ solving this
equation numerically. Figure 2 shows the behavior of the tangent plane for the
18-site honeycomb problem as a function of $\mu$ for select values of $U$ and
$\beta$. For large $\beta$ where the $\tanh$ becomes a sign function, the
tangent plane has plateaus that are connected by constant slopes. The location
of those depends on the spectrum of the hopping matrix $K$; an example is
shown in fig. 2.

####  II.4.1 Properties of the tangent plane in the $\mu=0$, $\beta\to\infty$
limit

The imaginary offset vanishes $\phi_{0}=0$ so that the tangent and real planes
coincide if the noninteracting energy eigenvalues $\epsilon_{k}$ are symmetric
about zero. This symmetry naturally occurs for bipartite lattices, such as the
honeycomb lattice, in the absence of chemical potential $\mu=0$. For these
cases $\phi_{0}=0$ for any inverse temperature $\beta$.

Moreover, in the $\beta\to\infty$ limit the $\tanh$ functions become the sign
function, and if there are equal numbers of positive and negative
noninteracting eigenvalues (not necessarily symmetric about zero), the sum
also vanishes and the tangent plane again corresponds to the real plane. Both
$C_{20}$ and $C_{60}$ have non-symmetric spectra but $C_{60}$ enjoys equally
many positive and negative noninteracting energies. In Figure 2 we show how
finite temperature smooths the piecewise-linear $\beta=\infty$ tangent plane
for the 18-site honeycomb lattice.

Figure 2: Tangent plane of 18-site honeycomb lattice. The vertical lines mark
transitions where the argument of a $\tanh$ that determines the tangent plane
(14) switches its sign. At the beginning of a downwards slope in
$\beta\to\infty$ the sign switches from $-1$ to $0$ and at the end from $0$ to
$+1$.

In fig. 3 we compare the behavior of the tangent plane for varying $U$, $\mu$,
$\beta$, and different lattices, both bipartite and non-bipartite. The results
at $\mu=0$ shown in this figure confirm our statements above. For the
fullerene results, which are non-bipartite, the choice of $\beta=10$ is large
enough that the resulting tangent plane at $\mu=0$ is nearly identical to the
real plane.

(a) Comparison of different bipartite lattices at $U=2$ and $\beta=10$. For
the bipartite lattices the tangent plane is antisymmetric around $\mu=0$.

(b) Comparison of the two fullerenes we investigated at $U=2$ and $\beta=10$.
For the non-bipartite lattices the tangent plane is not symmetric.

(c) The four site honeycomb lattice at $\beta=10$ for different $U$.

(d) The four site honeycomb lattice at $U=2$ for different $\beta$.

Figure 3: Tangent plane offsets (normalized to $\tilde{U}$) depending on the
chemical potential $\mu$, for various systems and parameters.

####  II.4.2 Properties of the tangent plane in the $\mu\beta\to\infty$ limit

Another interesting scenario is to consider the behavior of the tangent plane
in the $\mu\beta\to\infty$ limit. To understand the behavior in this limit we
start again with the action (9). With repeated application of Schur’s
complement, we can express (see ref. [9] for an explicit derivation)

$\displaystyle\log\det M[\pm]$
$\displaystyle=\log\det\left(\mathbbm{1}+\mathbb{F}[\pm]\right)$
$\displaystyle\mathbb{F}[\pm]=\prod_{t=0}^{N_{t}-1}[e^{\pm\tilde{K}}][e^{\mp
i\phi_{t}}]e^{\pm\tilde{\mu}}\ .$ (15)

where the 0th timeslice is rightmost and each term in $[$square brackets$]$ in
the product represents a space$\times$space matrix. Since the chemical
potential term is proportional to the identity, we can bring it out of the
product, so the action becomes

$S[\bm{\phi}]=\frac{\bm{\phi}^{2}}{2\tilde{U}}-\log\det\left(\mathbbm{1}+e^{-\mu\beta}\prod_{t=0}^{N_{t}-1}[e^{-\tilde{K}}][e^{i\phi_{t}}]\right)-\log\det\left(\mathbbm{1}+e^{+\mu\beta}\prod_{t=0}^{N_{t}-1}[e^{\tilde{K}}][e^{-i\phi_{t}}]\right)\
.$ (16)

In the limit of asymptotically large $\beta\mu$ the determinants simplify

$\displaystyle\lim_{\mu\beta\to\infty}S[\phi]=$
$\displaystyle\frac{\bm{\phi}^{2}}{2\tilde{U}}-\log\det\left(\mathbbm{1}\right)-\log\det\left(e^{+\mu\beta}\prod_{t=0}^{N_{t}-1}[e^{\tilde{K}}][e^{-i\phi_{t}}]\right)+\mathcal{O}(e^{-\mu\beta})$
$\displaystyle=$
$\displaystyle\frac{\bm{\phi}^{2}}{2\tilde{U}}-N_{x}\mu\beta-\log\det\left(\prod_{t=0}^{N_{t}-1}[e^{+\tilde{K}}][e^{-i\phi_{t}}]\right)$
$\displaystyle=$
$\displaystyle\frac{\bm{\phi}^{2}}{2\tilde{U}}-N_{x}\mu\beta-\log\left[e^{-i\Phi}\det
e^{+K\beta}\right]$ $\displaystyle=$
$\displaystyle\frac{\bm{\phi}^{2}}{2\tilde{U}}-N_{x}\mu\beta+i\Phi-\beta\tr{K}$
(17)

where we define $\Phi=\sum_{x,t}\phi_{x,t}$. Since our hopping matrices have
no self hopping $\tr{K}=0$. Now solving for critical points $\partial
S/\partial\phi=0$, we find the main critical point at constant field with
$\phi_{c}=-i\tilde{U}$. This demonstrates that the tangent plane approaches
$-\tilde{U}$ in the large $\mu\beta$ limit. Similarly, in the
$\mu\beta\to-\infty$ limit one finds

$S[\bm{\phi}]=\frac{\bm{\phi}^{2}}{2\tilde{U}}+N_{x}\mu\beta-i\Phi\ ,$ (18)

and the tangent plane approaches $+\tilde{U}$.

The resulting tangent plane shift in this limit has important implications for
our stochastic calculations. Simulating on the tangent plane means using
components of the field, $\phi_{j}$, that are offset by $-i\tilde{U}$, i.e.
$\phi_{j}\to\phi_{j}-i\tilde{U}\ \forall\ j$. The resulting action under this
deformation becomes

$S[\bm{\phi}-i\tilde{U}]=\frac{\bm{\phi}^{2}}{2\tilde{U}}-\frac{N_{x}U\beta}{2}+N_{x}U\beta-
N_{x}\mu\beta=\left(\frac{U}{2}-\mu\right)\beta
N_{x}+\frac{\bm{\phi}^{2}}{2\tilde{U}}\ .$ (19)

This means that the action on this flat contour is purely real in this limit.
That is, the tangent plane _completely solves the sign problem_ in this limit.
It is equivalent, up to some overall shift in the energy, to a _quenched_
calculation, with no fermion matrix. Later we show numerical results that
confirm these findings.

###  II.5 Quantum corrections to the saddle point

The saddle point that defines the tangent plane corresponds to the critical
point of the classical action. In quantum field theory the location of this
point shifts due to the presence of quantum fluctuations which, in our case,
corresponds to thermal fluctuations. We can estimate this shift by calculating
the quantum effective action and determining the extremum of this action, as
is done in standard textbooks on QFT. This correction to the saddle point
corresponds to the inclusion of all one-particle irreducible (1PI) diagrams.
Thus it represents a quantum (thermal) correction to the classical saddle
point, and the ensuing constant manifold that intercepts this point is
expected to reduce the sign problem. We assume that the maximum we find, when
including higher order correction terms, will return an offset that reduces
the sign-problem even more than the basic tangent plane. To start, we assume
that $\phi_{c}$ is the saddle point in the presence of quantum fluctuations
and apply the saddle point approximation about this point. That is, we expand
the action in powers of a small perturbation $\eta$ about this point,
$\phi=\phi_{c}+\eta$. This gives

$\displaystyle S[\phi_{c}+\bm{\eta}]=$ $\displaystyle
S[\phi_{c}]+\left(\bm{\eta}\cdot\bm{\nabla}\right)S[\phi_{c}]+\frac{1}{2}\left(\bm{\eta}\cdot\bm{\nabla}\right)^{2}S[\phi_{c}]+\order{\eta^{3}}$
$\displaystyle=$ $\displaystyle
S[\phi_{c}]+\left(\bm{\eta}\cdot\bm{\nabla}\right)S[\phi_{c}]+\frac{1}{2}\bm{\eta}\cdot\mathbb{H}_{S[\phi_{c}]}\cdot\bm{\eta}+\order{\eta^{3}}.$
(20)

Here we have made use of the Hessian,

$\left(\mathbb{H}_{S[\phi_{c}]}\right)_{x^{\prime}t^{\prime},xt}=\left(\partial_{x^{\prime}t^{\prime}}\partial_{xt}S[\phi]\right)\vline_{\phi=\phi_{c}}\
.$ (21)

Since $\eta$ is assumed small, we will omit the $\order{\eta^{3}}$ terms.
Furthermore because the critical point satisfies $\evaluated{\nabla
S[\phi]}_{\phi_{c}}=0$ the linear terms also vanish. Hence the path integral
simplifies to a Gaussian integral which we can do,

$\int\mathcal{D}\phi\,e^{-S[\phi]}\approx
e^{-S[\phi_{c}]}\int\mathcal{D}\bm{\eta}\,e^{-\frac{1}{2}\bm{\eta}\cdot\mathbb{H}_{S[\phi_{c}]}\cdot\bm{\eta}}=e^{-S[\phi_{c}]}\left(\det\mathbb{H}_{S[\phi_{c}]}\right)^{-1/2}\equiv
e^{-S_{\rm eff}[\phi_{c}]}\ .$ (22)

This allows us to formulate an effective action

$S_{\rm
eff}[\phi_{c}]=S[\phi_{c}]+\frac{1}{2}\log\det\mathbb{H}_{S[\phi_{c}]}\ .$
(23)

The extremum of this action defines our 1PI-corrected spacetime-constant
saddle point $\phi_{{c}}=i\phi_{1}$. Note that without the Hessian term we
recover our original action and the extremum in this case is the saddle point
of our leading order classical action that defines the tangent plane (14). In
comparison, our 1PI-corrected effective action (23) includes the quantum
effects at next to leading order (NLO). In Appendix C we show how to evaluate
the Hessian (21) when $\phi_{{c}}=i\phi_{1}$ a spacetime constant.

###  II.6 Excursion to infinite lattices

In this section we demonstrate on select lattices how to determine the tangent
plane in the infinite-volume limit. We provide two well known examples: the
2-dimensional square and honeycomb lattices. In the infinite volume limit we
can access every mode in the first Brillouin zone (B.Z.) and we can replace
the sum over noninteracting energies in eq. 14 with a momentum integral. For a
2-dimensional square lattice one has

$\frac{1}{N_{x}}\sum_{k}\to\int_{\bm{k}\in B.Z.}\frac{d\bm{k}}{(2\pi)^{2}},$

where $\bm{k}\equiv(k_{x},k_{y})$ with $-\pi\leq k_{i}<\pi$ (square B.Z.). The
non-interacting energies are given by

$\epsilon(\bm{k})=2\left(\cos(k_{x})+\cos(k_{y})\right)\ .$ (24)

Making these substitutions to determine the tangent plane (14) leads to

$\phi_{0}/\delta=-U\int_{\bm{k}\in
B.Z.}\frac{d\bm{k}}{(2\pi)^{2}}\tanh\left(\frac{1}{2}\beta\left[\epsilon(\bm{k})+\mu+\phi_{0}/\delta\right]\right)\
.$ (25)

For the infinite honeycomb lattice the non-orthogonal lattice translation
vectors and the two-band structure means we must substitute

$\frac{1}{N_{x}}\sum_{k}\to\frac{3\sqrt{3}}{2}\int_{\bm{k}\in
B.Z.}\frac{d\bm{k}}{(2\pi)^{2}}\frac{1}{2}\sum_{\sigma=\pm 1}\ ,$

where the factor $3\sqrt{3}/2$ comes from the hexagonal geometry of the B.Z.
and $\sigma$ runs over the two bands. The non-interacting energies are [35,
36]

$\displaystyle\epsilon_{\pm}(\bm{k})$ $\displaystyle=\pm|f(\bm{k})|$
$\displaystyle f(\bm{k})$
$\displaystyle=1+2e^{-\frac{3ik_{x}}{2}}\cos\left(\frac{\sqrt{3}k_{y}}{2}\right)\
.$ (26)

So we find

$\phi_{0}/\delta=-U\frac{3\sqrt{3}}{2}\int_{\bm{k}\in
B.Z.}\frac{d\bm{k}}{(2\pi)^{2}}\frac{1}{2}\sum_{\sigma=\pm
1}\tanh\left(\frac{1}{2}\beta\left[\epsilon_{\sigma}(\bm{k})+\mu+\phi_{0}/\delta\right]\right)\
.$ (27)

Figure 4: The zero-temperature tangent plane $\phi_{0}$ (normalized to
$\tilde{U}$) for the infinite 2-D square lattice (left) and honeycomb lattice
(right) as a function of chemical potential $\mu$ for various onsite
interactions $U$, as labelled in the figure.

We solve the square-lattice (25) and honeycomb (27) relations numerically. In
fig. 4 we show the solutions of $\phi_{0}$ for both the square and honeycomb
system for select values of $U$. We see that $\phi_{0}$ remains smooth as a
function of chemical potential $\mu$, even in the limit of zero temperature
($\beta\gg 1$). Also, in both cases, in the limit of asymptotically large
$\mu$ we have $\phi_{0}\to-\tilde{U}$.

##  III Numerical Optimization Method

In many cases both tangent plane and NLO offsets lead to an improvement in
statistical power, with NLO typically providing modest improvement over the
tangent plane (but not in all cases). However, a simple numerical
investigation shows that one can further improve the statistical power, in
most cases by shifting beyond the NLO result. For example, in fig. 5 we show
the statistical power for the eight-site honeycomb system coming from a scan
of various offsets that include the real plane, tangent plane, and the NLO
offset. The scan shows a singular peak in the statistical power. However, this
peak does not occur at either the tangent plane or the NLO offset.

Figure 5: Visualization of effect of an imaginary offset on the sign problem.
This example shows the eight site honeycomb lattice with $N_{t}=16$,
$\beta=8$, $U=2$ and $\mu=1$. The errorbands show the bootstrap errors of the
first and second derivative at the datapoint.

In all our investigations of different systems to date we find similar
behavior; namely, there is a singular peak in statistical power due to
constant offset. We refer to this offset that maximizes the statistical power
as the _optimized shift_.

Because the greater statistical power means smaller required samples (8), the
potential savings in computational resources when simulating at the optimized
shift in comparison to either the NLO or tangent plane can be orders of
magnitude. However, determining the location of the peak from a simple raster
scan in offsets is timely and inefficient, since each point requires an HMC
simulation with sufficient statistics to resolve the statistical power.
Instead, we formulate a search algorithm like Newton-Raphson, relying on the
calculation of derivatives of the statistical power using the _current_ HMC
ensemble to make a prediction for the location of the offset that corresponds
to the peak of the statistical power. We then iterate this procedure to
converge to the peak.

Our algorithm requires the first two derivatives of the statistical power with
respect to the imaginary offset $\phi_{0}$ from the existing Markov Chain
configurations,

$\displaystyle\derivative{\phi_{0}}\expectationvalue{e^{-iS_{I,\phi_{0}}}}_{R,\phi_{0}}=$
$\displaystyle\derivative{\phi_{0}}\frac{\int\mathcal{D}\phi
e^{-S_{\phi_{0}}}}{\int\mathcal{D}\phi e^{-S_{R,{\phi_{0}}}}}$
$\displaystyle=$
$\displaystyle\expectationvalue{e^{-iS_{I,{\phi_{0}}}}}_{R,{\phi_{0}}}\expectationvalue{\derivative{S_{R,{\phi_{0}}}}{{\phi_{0}}}}_{R,{\phi_{0}}}-\expectationvalue{e^{-S_{I,{\phi_{0}}}}\derivative{S_{{\phi_{0}}}}{{\phi_{0}}}}_{R,{\phi_{0}}}$
(28)
$\displaystyle\derivative{{}^{2}}{{\phi_{0}}^{2}}\expectationvalue{e^{-iS_{I,{\phi_{0}}}}}_{R,{\phi_{0}}}=$
$\displaystyle\derivative{\phi_{0}}\left(\expectationvalue{e^{-iS_{I,{\phi_{0}}}}}_{R,{\phi_{0}}}\expectationvalue{\derivative{S_{R,{\phi_{0}}}}{{\phi_{0}}}}_{R,{\phi_{0}}}-\expectationvalue{e^{-S_{I,{\phi_{0}}}}\derivative{S_{{\phi_{0}}}}{{\phi_{0}}}}_{R,{\phi_{0}}}\right)$
$\displaystyle=$
$\displaystyle\expectationvalue{e^{-iS_{I,{\phi_{0}}}}}_{R,{\phi_{0}}}\left(2\expectationvalue{\derivative{S_{R,{\phi_{0}}}}{{\phi_{0}}}}_{R,{\phi_{0}}}^{2}+\expectationvalue{\derivative[2]{S_{R,{\phi_{0}}}}{{\phi_{0}}}-\derivative{S_{R,{\phi_{0}}}}{{\phi_{0}}}^{2}}_{R,{\phi_{0}}}\right)$
$\displaystyle-2\expectationvalue{e^{-S_{I,{\phi_{0}}}}\derivative{S_{{\phi_{0}}}}{{\phi_{0}}}}_{R,{\phi_{0}}}\expectationvalue{\derivative{S_{R,{\phi_{0}}}}{{\phi_{0}}}}_{R,{\phi_{0}}}-\expectationvalue{e^{-S_{I,{\phi_{0}}}}\left(\derivative[2]{S_{{\phi_{0}}}}{{\phi_{0}}}-\derivative{S_{{\phi_{0}}}}{{\phi_{0}}}^{2}\right)}_{R,{\phi_{0}}}$
(29)

We stress that the calculations of these derivatives rely only on a single
ensemble. Ref. [37] points out that these derivatives may be simplified and
estimated with reliability even in cases with a sign problem. Practically,
they enable iterative procedures for a predicting the optimized shift. When
the second derivative is negative we can get a good prediction via Newton-
Raphson

$\phi_{0,i+1}=\phi_{0,i}-\frac{\derivative{\phi_{0}}\expectationvalue{e^{-iS_{I,\phi_{0}}}}_{R,\phi_{0},i}}{\derivative{{}^{2}}{{\phi_{0}}^{2}}\expectationvalue{e^{-iS_{I,{\phi_{0}}}}}_{R,{\phi_{0}},i}}\
.$ (30)

When the second derivative is positive we work with the first derivative to
approach the peak. As soon as there are points with opposite first
derivatives, their central value usually gets us into the region where the
second derivative is negative or at least fairly close to it. Additionally we
limit the searching region to the interval $[-\tilde{U},+\tilde{U}]$. In
principle, higher order derivatives can also be calculated and used to predict
the location of the optimized shift though the statistical errors in these
terms grow. In practice we find that the procedure converges quickly when
starting from a region where the sign problem is light enough to calculate at
least the first derivative. However it can fail when this is not the case and
it gets stuck when the statistical power is of the same order of magnitude as
its uncertainty. A more quantitative demonstration of this method will be
given in section III.1. A more advanced method might fit all known
measurements and derivatives to estimate the location of the peak.

When we are only interested in separate sets of parameters we have to rely on
an analytic approximation as the starting point. When we want to scan over one
parameter sufficiently finely, we can do so iteratively starting from the
previous offset or a rescaled version of it. In this paper we rescaled it with
the fraction of new and previous tangent plane offset.

Figure 6 gives a cartoon showing the different offsets in comparison to the
Lefschetz thimbles and conveys a geometrically intuitive understanding as to
why some planes do better than others.

$\Im{\phi}$$\Re{\phi}$Lefschetz ThimbleReal PlaneTangent PlaneNLO
PlaneOptimized Plane

Figure 6: A cartoon of the different manifolds referred to throughout this
work. The Lefschetz thimbles are drawn to resemble contours of holomorphic
flow applied to constant fields. The real plane, tangent plane, NLO estimate
(next to leading order correction) and optimized plane show the planar
manifolds that we use as our integration regions. The red dot marks the main
critical point.  Figure 7: The datapoints in this heatmap show the statistical
power of the system as indicated by the colorbar at different chemical
potential $\mu$ evaluated with HMC on a plane with an imaginary offset given
by the y-axis. The red and purple curve show our analytically determined
offsets. The blue curve connects the offsets with the greatest statistical
power at each $\mu$. The inset plot shows a slice of the heatmap to visualize
the connection to fig. 5. The system is the eight site honeycomb lattice with
$N_{t}=16$, $\beta=8$, $U=2$.

In fig. 7 we show the statistical power for an interesting range of chemical
potentials and imaginary offsets. We trace contours for the tangent plane, NLO
offset, and the best offsets. The best case scenario would be a cheaply
determined estimate of the optimized offset for a given $\mu$. The statistical
powers of the real-, tangent- , NLO, and optimized planes are compared as a
function of $\mu$ in fig. 8. The key takeaways are that the tangent plane
consistently and drastically outperforms the standard algorithm at practically
the same cost, and that the sign problem vanishes when the system becomes
saturated. Furthermore we see that we can reduce the sign problem for
parameters where the tangent plane is insufficient. For system sizes of
interest these differences determine whether a system can be calculated (with
reasonable resources) or not.

Figure 8: Comparing statistical power of the real plane, tangent plane, NLO
correction and optimized shift as a function of $\mu$. The system is an eight
site honeycomb lattice with $N_{t}=16$, $\beta=8$ and $U=2$.

Note that the sign problem also gets worse with increasing $\beta$ and $U$
such that even the optimized plane will eventually fail at finite $\mu$. This
simple method is not suitable to fix the sign problem across the board, it
just leads to an efficient expansion of the calculable parameter space, which
might or might not include interesting physical phenomena, but definitely
enables us to do better zero temperature extrapolations. For stronger sign
problems we have to rely on manifolds with more parameters, either with simple
parametrizations or with neural networks [15, 14, 13].

###  III.1 Benchmarks

In this section we demonstrate the convergence of our numerical optimization
algorithm and present the resulting increase of the statistical power. The
examples refer to the eight site honeycomb lattice with $N_{t}=16$, $\beta=8$
and $U=2$. Figure 9 showcases the convergence with iterations from different
starting points. While they would all converge to the same offsets eventually,
we observe that starting in a region with unclear first derivative, that is
dominated by statistical noise, turns the algorithm into a random walk, which
can be observed in the real plane example. This further highlights the
importance of having good analytic starting points.

(a) Starting from real plane.

(b) Starting from tangent plane.

(c) Starting from leading order correction.

(d) Iteratively starting at rescaled previous result.

Figure 9: Numerically optimized offset after a given number of iterations.
Iteration 0 marks the starting offsets. This plot shows the convergence to the
optimal offset and the value of a good starting guess. The iterative method
started from $\mu=0$ and $\mu=6$, meeting in the middle.

Figure 10 shows the significant improvements in statistical power that can be
achieved by just a few iterations. We observe in both figures that most runs
starting from a reasonable guess converge and roughly agree with each other
after just 3 iterations. Here the improvements from tangent plane to leading
order correction to iterative starting points can be best observed by
comparing the second iterations with each other and and in comparison to the
converged result.

(a) Starting from real plane.

(b) Starting from tangent plane.

(c) Starting from leading order correction.

(d) Iteratively starting at rescaled previous result.

Figure 10: Statistical power at numerically determined optimal offset, after a
given number of optimization steps. For each iteration we used $10\,000$ HMC
steps (+tuning).

##  IV Results

In this section we provide physical observables determined by HMC with the
introduced modifications. The observables of our choice are the single
particle correlation functions

$C_{k}(\tau)=\expectationvalue{a_{k}(\tau)a_{k}^{\dagger}(0)}$ (31)

where $a$ and $a{}^{\dagger}$ are particle ladder operators and $k$ labels
operators in definite the irreducible representations of the lattice
automorphism group. In the honeycomb case $k$ labels operators with definite
momentum; for the fullerenes $k$ labels representations of the icosahedral
symmetry group. We also sum the local charges (1) to measure the global charge

$\expectationvalue{Q}=\expectationvalue{\sum_{x}q_{x}}=N_{x}-2\sum_{k}C_{k}(\tau=0)\
.$ (32)

To establish trust in the physical correctness of the algorithm we compare
with the eight site honeycomb lattice for which we have exact results from
direct diagonalization. We find that the uncertainties of the standard real
plane HMC are much greater than the uncertainties of our more advanced
methods. Especially for the eight site lattice the real plane results do not
agree with the exact solution while the calculations with an imaginary offset
match it very well. Above all the optimized offset resembles the exact
solution with great precision. This plus the agreement with the NLO, and often
with the tangent plane as well, gives us confidence in the results of the
larger systems. We present the correlation functions resulting from our
methods in fig. 11.

(a) Eight site honeycomb lattice at $\beta=8$.

(b) 18 site honeycomb lattice at $\beta=8$.

(c) $C_{20}$ at $\beta=6$.

(d) $C_{60}$ at $\beta=6$.

Figure 11: Single particle correlation function for different lattices and
calculated with different methods. All systems were evaluated with $N_{t}=32$
at $U=2$ and $\mu=1$. Each data point was calculated from a markov chain with
$50,000$ HMC configurations, where we measured on each $10^{\text{th}}$ to
reduce autocorrelation, furthermore we averaged over correlators guaranteed to
be equal by symmetry. The exact solution in fig. 11(a) was determined by exact
diagonalization in the temporal continuum limit.

Figure 12 shows that for the same number of configurations the quality of the
measured observables seems to match the expected outcome from comparing the
statistical powers, which can be found in fig. 8 and fig. 13. We see larger
statistical fluctuations with worse statistical power; the optimized method
consistently performs best. Figure 12(a) shows that most of the numerical
results estimate the charge correctly (32) according to the exact results. The
real plane HMC underestimates its error systematically for large ranges of
$\mu$. Still, the optimized offset resembles the exact result best; the
tangent plane and NLO arguably offer comparable uncertainties for many
parameter choices. Figure 12(d) and fig. 13(c) show that our method has limits
and not every sign problem can be conquered with a simple constant offset.
Also in certain areas of the other charge plots, we see that the optimization
routine could fail when the sign problem is very strong, causing a worse
result than the NLO. An interesting observation is that where the charge
flattens the statistical power of the tangent and NLO planes peak.

(a) Eight site honeycomb lattice at $\beta=8$.

(b) 18 site honeycomb lattice at $\beta=8$.

(c) $C_{20}$ at $\beta=6$.

(d) $C_{60}$ at $\beta=6$.

Figure 12: Charge expectation value at varying $\mu$ for different lattices
and calculated with different methods. All systems have $N_{t}=16$ and $U=2$.
Each data point was calculated from a markov chain with $50\,000$ HMC
configurations, where we measured on each $10^{\text{th}}$ to reduce
autocorrelation. The exact solution in fig. 12(a) was determined by exact
diagonalization with $N_{t}=16$ discretization.

##  V Conclusion

Our results clearly show that the simple introduction of an imaginary shift of
the integration contour can greatly impact the severity of the sign problem in
quantum field theory with practically no additional computation cost or human
effort required. We provide two analytic expressions for such offsets for the
Hubbard model. Further a careful tuning of these offsets while coming at a
small cost can lead to even better outcomes, especially when a range of
parameters is to be explored. The reduction of the sign problem depends on the
system and the physical parameters, but can potentially make a difference of
orders of magnitude which can be seen as an increase in measurement precision
at a fixed sample size or saving of computational resources on the way to
achieve a certain desired precision. Even though this does not eliminate the
sign problem entirely, it extends the parameter space that is explorable
within our naturally limited resources and allows us to perform higher quality
extrapolations.

In addition to the analysis of the method itself we provided observables to
condensed matter systems which could not be calculated before with lattice
stochastic methods, going as far as the $C_{60}$ buckyball a stable
synthesizable carbon nanosystem. Further we found that in the limit of
infinite chemical potential, which in case of the Hubbard model refers to a
completely filled/empty lattice, the sign problem vanishes at a certain offset
stating the question if this behavior could be seen in other theories as well.
Our numerical optimization could lend itself to unsupervised learning, driving
towards maximizing the statistical power and minimizing its derivative (III).

There are some open questions remaining in regards to combining optimized
offset with neural networks, that we will address in the future. Is the
optimized offset also the best starting point for generating training data?
Does the uplift of neural networks and optimized shift over the tangent plane
correlate? Further we plan to investigate the $C_{60}$ lattice in more detail,
as our methods have opened the door to high quality measurements on these
large systems that are out of reach for exact diagonalization as well as HMC
without sign problem optimization. We are developing a library for these
nanosystems calculations with the intention of making it publicly available so
everyone can try out our methods on the theory or model they are interested
in.

###### Acknowledgements.

We thank Neill Warrington for many helpful discussions related to this work as
well as Timo Lähde for his valuable comments. This work was funded in part by
the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through
the funds provided to the Sino-German Collaborative Research Center
“Symmetries and the Emergence of Structure in QCD” (NSFC Grant No.
12070131001, DFG Project-ID 196253076 – TRR110) as well as the STFC
Consolidated Grant ST/T000988/1. We gratefully acknowledge the computing time
granted by the JARA Vergabegremium and provided on the JARA Partition part of
the supercomputer JURECA at Forschungszentrum Jülich.

##  Appendix A Statistical Power

In fig. 13 we show the behavior of the statistical power over the interesting
range of $\mu$ for the remaining lattices, where we again compare the
different offset methods. We also see in the plot that for non-bipartite
lattices introducing a small chemical potential can relieve the sign problem.

(a) 18 site honeycomb lattice with $\beta=8$.

(b) $C_{20}$ lattice with $\beta=6$.

(c) $C_{60}$ lattice with $\beta=6$.

Figure 13: Statistical power for different lattices over a range of $\mu$. For
all systems $U=2$ and $N_{t}=16$.

##  Appendix B Detailed derivation of tangent plane

Starting from the action (9) we provide a detailed derivation the equation
that determines the tangent plane (14). As in the main text (15) let

$\displaystyle\mathbb{F}_{\pm}=\mathbb{F}[\pm\phi,\pm K,\pm\mu]$
$\displaystyle=\prod_{\tau=0}^{N_{t}-1}f^{\pm}_{\tau}$ $\displaystyle
f^{\pm}_{\tau}=\left[e^{\pm\tilde{K}}\right]\left[e^{\pm(-i\phi_{x,\tau}+\tilde{\mu})}\right]$
(33)

where each term in $[$square brackets$]$ is a space-by-space matrix and the
product is from right ($\tau=0$) to left $(\tau=N_{t}-1)$. Then, using the
Schur complement

$\displaystyle\det
M[\phi,K,\mu]=\det\left(\mathbbm{1}+\mathbb{F}[\phi,K,\mu]\right)$ (34)

and likewise for $M[-\phi,-K,-\mu]$. Since the inverses exist for any
configuration with finite action

$\displaystyle\frac{\partial S}{\partial\phi_{x,t}}$
$\displaystyle=\frac{1}{\tilde{U}}\phi_{x,t}-\tr{(\mathbbm{1}+\mathbb{F}_{+})^{-1}\partial_{\phi_{x,t}}\mathbb{F}_{+}}-\tr{(\mathbbm{1}+\mathbb{F}_{-})^{-1}\partial_{\phi_{x,t}}\mathbb{F}_{-}}$
(35)

where the two traces correspond to the particle and hole fermion matrices
(15).

Each auxiliary field only appears once in any given $\mathbb{F}$, inside one
entry of a diagonal matrix. So, differentiating $\mathbb{F}_{\pm}$ inserts
$\mp i\mathbb{P}_{x}$ where $\mathbb{P}_{x}$ projects to site $x$,

$\displaystyle\partial_{\phi_{x,t}}\mathbb{F}_{\pm}$
$\displaystyle=\left[\prod_{\tau=t}^{N_{t}-1}f^{\pm}_{\tau}\right]\left[\mp
i\mathbb{P}_{x}\right]\left[\prod_{\tau=0}^{t-1}f^{\pm}_{\tau}\right]$ (36)

where the products go from right to left. Cycling the traces so that the
projector is rightmost gives

$\displaystyle\frac{\partial S}{\partial\phi_{x,t}}$
$\displaystyle=\frac{1}{\tilde{U}}\phi_{x,t}+i\sum_{s}s\tr{\left[\prod_{\tau=0}^{t-1}f^{\pm}_{\tau}\right](\mathbbm{1}+\mathbb{F}_{s})^{-1}\left[\prod_{\tau=t}^{N_{t}-1}f^{\pm}_{\tau}\right]\mathbb{P}_{x}}$
(37)

where the sum over $s$ runs over $\pm 1$.

Plugging in the imaginary spacetime constant $\phi=i\phi_{0}$ means the
auxiliary field factors are proportional to the identity matrix,

$\displaystyle
f^{\pm}_{\tau}=\exp\left(\pm(\tilde{K}+\phi_{0}+\tilde{\mu})\right)$ (38)

independent of timeslice $\tau$. Since $\phi_{0}$ is space-independent we can
sum on $x$ and use the completeness $\sum_{x}\mathbb{P}_{x}=\mathbbm{1}$.
Since $N_{t}\tilde{K}=\beta K$ and $N_{t}\tilde{\mu}=\beta\mu$,

$\displaystyle 0=\frac{1}{\tilde{U}}N_{x}\phi_{0}+\sum_{s\in{\pm
1}}s\tr{\frac{e^{s(\beta K+N_{t}\phi_{0}+\beta\mu)}}{1+e^{s(\beta
K+N_{t}\phi_{0}+\beta\mu)}}}$ (39)

and evaluating this relationship in the eigenbasis of $K$ yields the tangent
plane relation (14).

##  Appendix C Detailed derivation of NLO

To find the NLO constant imaginary offset we need to minimize the effective
action (23) which requires computing the Hessian

$\displaystyle\mathbb{H}_{x^{\prime}t^{\prime},xt}=\left(\partial_{x^{\prime}t^{\prime}}\partial_{xt}S[\phi]\right)\vline_{\phi=i\phi_{1}}\
.$ (40)

We start from the general single derivative (35) and differentiate again.
Without loss of generality we assume $t^{\prime}\geq t$,

$\displaystyle\partial_{x^{\prime}t^{\prime}}\partial_{xt}S[\phi]=\frac{1}{\tilde{U}}\delta_{x^{\prime}x}\delta_{t^{\prime}t}$
$\displaystyle-\sum_{s}\tr{(\mathbbm{1}+\mathbb{F}_{s})^{-1}\partial_{\phi_{x^{\prime},t^{\prime}}}\partial_{\phi_{x,t}}\mathbb{F}_{s}-(\mathbbm{1}+\mathbb{F}_{s})^{-1}\left(\partial_{\phi_{x^{\prime},t^{\prime}}}\mathbb{F}_{s}\right)(\mathbbm{1}+\mathbb{F}_{s})^{-1}\left(\partial_{\phi_{x,t}}\mathbb{F}_{s}\right)}$
(41)

The second derivative of $\mathbb{F}$ is much like the first (36) but with a
second projector $\mathbb{P}_{x^{\prime}}$ inserted at time $t^{\prime}$; the
$(\mp i)^{2}=-1$ regardless of sign choice.

$\displaystyle\partial_{x^{\prime}t^{\prime}}\partial_{xt}\mathbb{F}_{s}$
$\displaystyle=-\left[\prod_{\tau=t^{\prime}}^{N_{t}-1}f^{s}_{\tau}\right]\mathbb{P}_{x^{\prime}}\left[\prod_{\tau=t}^{t^{\prime}-1}f^{s}_{\tau}\right]\mathbb{P}_{x}\left[\prod_{\tau=0}^{t-1}f^{s}_{\tau}\right].$
(42)

In fact, this looks much like the two-inverse term, though that term has an
inverse between the projectors,

$\displaystyle\left(\partial_{\phi_{x^{\prime},t^{\prime}}}\mathbb{F}_{s}\right)(\mathbbm{1}+\mathbb{F}_{s})^{-1}\left(\partial_{\phi_{x,t}}\mathbb{F}_{s}\right)$
$\displaystyle=-\left[\prod_{\tau=t^{\prime}}^{N_{t}-1}f^{s}_{\tau}\right]\mathbb{P}_{x^{\prime}}\left[\prod_{\tau=0}^{t^{\prime}-1}f^{s}_{\tau}\right]\left[\mathbbm{1}+\mathbb{F}_{s}\right]^{-1}\left[\prod_{\tau=t}^{N_{t}-1}f^{s}_{\tau}\right]\mathbb{P}_{x}\left[\prod_{\tau=0}^{t-1}f^{s}_{\tau}\right].$
(43)

Cycling the trace so that $\mathbb{P}_{x}$ is rightmost and consolidating like
factors gives

$\displaystyle\partial_{x^{\prime}t^{\prime}}\partial_{xt}S[\phi]=\frac{1}{\tilde{U}}\delta_{x^{\prime}x}\delta_{t^{\prime}t}$
$\displaystyle+\sum_{s}\tr{\left[\prod_{\tau=0}^{t-1}f^{s}_{\tau}\right](\mathbbm{1}+\mathbb{F}_{s})^{-1}\left[\prod_{\tau=t^{\prime}}^{N_{t}-1}f^{s}_{\tau}\right]\mathbb{P}_{x^{\prime}}\left[\prod_{\tau=t}^{t^{\prime}-1}f^{s}_{\tau}\right]\left[\mathbbm{1}-\left[\prod_{\tau=0}^{t-1}f^{s}_{\tau}\right]\left[\mathbbm{1}+\mathbb{F}_{s}\right]^{-1}\left[\prod_{\tau=t}^{N_{t}-1}f^{s}_{\tau}\right]\right]\mathbb{P}_{x}}.$
(44)

Making repeated use of $C^{-1}B^{-1}A^{-1}=(ABC)^{-1}$ we can re-express the
term in the sum

$\displaystyle\left[\prod_{\tau=0}^{t-1}f^{s}_{\tau}\right]\left[\mathbbm{1}+\mathbb{F}_{s}\right]^{-1}\left[\prod_{\tau=t}^{N_{t}-1}f^{s}_{\tau}\right]$
$\displaystyle=\left[\prod_{\tau=t-1}^{0}(f^{s}_{\tau})^{-1}\right]^{-1}\left[\mathbbm{1}+\prod_{\tau=0}^{N_{t}-1}f^{s}_{\tau}\right]^{-1}\left[\prod_{\tau=N_{t}-1}^{t}(f^{s}_{\tau})^{-1}\right]^{-1}$
$\displaystyle=\left[\left[\prod_{\tau=N_{t}-1}^{t}(f^{s}_{\tau})^{-1}\right]\left[\mathbbm{1}+\prod_{\tau=0}^{N_{t}-1}f^{s}_{\tau}\right]\left[\prod_{\tau=t-1}^{0}(f^{s}_{\tau})^{-1}\right]\right]^{-1}$
$\displaystyle=\left[\left[\prod_{\tau=t-1}^{t}\left(f^{s}_{\tau}\right)^{-1}\right]+\mathbbm{1}\right]^{-1}$
(45)

where the products of $f^{-1}$s count down from right to left and wrap from
$0$ to $N_{t}-1$. Inserting convenient products equal to the identity to use
the result (45) twice gives

$\displaystyle\partial_{x^{\prime}t^{\prime}}\partial_{xt}S[\phi]=$
$\displaystyle\frac{1}{\tilde{U}}\delta_{x^{\prime}x}\delta_{t^{\prime}t}$
(46)
$\displaystyle+\sum_{s}\tr{\left[\mathbbm{1}+\prod_{\tau=t-1}^{t}\left(f^{s}_{\tau}\right)^{-1}\right]^{-1}\left[\prod_{\tau=t^{\prime}-1}^{t}(f^{s}_{\tau})^{-1}\right]\mathbb{P}_{x^{\prime}}\left[\prod_{\tau=t}^{t^{\prime}-1}f^{s}_{\tau}\right]\left[\prod_{\tau=t-1}^{t}\left(f^{s}_{\tau}\right)^{-1}\right]\left[\mathbbm{1}+\prod_{\tau=t-1}^{t}\left(f^{s}_{\tau}\right)^{-1}\right]^{-1}\mathbb{P}_{x}}.$

Casting the inverse factors after $\mathbb{P}_{x^{\prime}}$ into the
denominator we arrive at

$\displaystyle\partial_{x^{\prime}t^{\prime}}\partial_{xt}S[\phi]=\frac{1}{\tilde{U}}\delta_{x^{\prime}x}\delta_{t^{\prime}t}$
$\displaystyle+\sum_{s}\tr{\left[\mathbbm{1}+\prod_{\tau=t-1}^{t}\left(f^{s}_{\tau}\right)^{-1}\right]^{-1}\left[\prod_{\tau=t^{\prime}-1}^{t}(f^{s}_{\tau})^{-1}\right]\mathbb{P}_{x^{\prime}}\left[\prod_{\tau=t}^{t^{\prime}-1}f^{s}_{\tau}\right]\left[\mathbbm{1}+\prod_{\tau=t}^{t-1}f^{s}_{\tau}\right]^{-1}\mathbb{P}_{x}},$
(47)

a convenient general form for arbitrary $\phi$.

To evaluate the effective action (23) and find the NLO imaginary offset we set
$\phi=i\phi_{1}$ a spacetime constant. Like before (38), the auxiliary field
terms become proprotional to the identity matrix and we can group terms into
powers of

$\displaystyle f^{\pm}=f^{\pm}_{\tau}=\exp\pm(\delta K+\delta\mu+\phi_{1})$
(48)

which has the nice property $(f^{\pm})^{-1}=f^{\mp}$ so that we may treat the
sign label as a true exponent. Defining $\Delta t=t^{\prime}-t$ we simplify to

$\displaystyle\mathbb{H}_{x^{\prime}t^{\prime},xt}=\left.\partial_{x^{\prime}t^{\prime}}\partial_{xt}S[\phi]\right|_{\phi=i\phi_{1}}=\frac{1}{\tilde{U}}\delta_{x^{\prime}x}\delta_{t^{\prime}t}$
$\displaystyle+\sum_{s}\tr{\left[\mathbbm{1}+f^{-sN_{t}}\right]^{-1}f^{-s\Delta
t}\mathbb{P}_{x^{\prime}}f^{+s\Delta
t}\left[\mathbbm{1}+f^{+sN_{t}}\right]^{-1}\mathbb{P}_{x}}.$ (49)

Since the hopping amplitudes are symmetric $K=K{}^{\top}$, the matrices $f$
(48) are too $f=f{}^{\top}$. Moreover, the projectors are symmetric
$\mathbb{P}=\mathbb{P}{}^{\top}$. Because the sum is over $s\in\\{\pm 1\\}$,
the signs and inverses conspire so that the two traces are over a matrix and
its transpose, and are therefore equal. So, we can consolidate the traces and
use the projectors to isolate needed matrix elements

$\displaystyle\mathbb{H}_{x^{\prime}t^{\prime},xt}=\frac{1}{\tilde{U}}\delta_{x^{\prime}x}\delta_{t^{\prime}t}$
$\displaystyle+2\tr{(\mathbbm{1}+f^{-N_{t}})^{-1}f^{-\Delta
t}\mathbb{P}_{x^{\prime}}f^{+\Delta
t}(\mathbbm{1}+f^{+N_{t}})^{-1}\mathbb{P}_{x}}$ (50)
$\displaystyle=\frac{1}{\tilde{U}}\delta_{x^{\prime}x}\delta_{t^{\prime}t}$
$\displaystyle+2\left[(\mathbbm{1}+f^{-N_{t}})^{-1}f^{-\Delta
t}\right]_{xx^{\prime}}\left[f^{+\Delta
t}(\mathbbm{1}+f^{+N_{t}})^{-1}\right]_{x^{\prime}x}.$ (51)

We may quickly evaluate the matrix elements by a unitary transformation from
the eigenbasis of $K$,

$\displaystyle\left[f^{\pm\Delta t}(\mathbbm{1}+f^{\pm
N_{t}})^{-1}\right]_{x^{\prime}x}$
$\displaystyle=\sum_{k}\mathbb{U}{}^{\dagger}_{x^{\prime}k}\frac{e^{\pm\Delta
t(\delta\epsilon_{k}+\phi_{1}+\delta\mu)}}{1+e^{\pm\beta(\epsilon_{k}+\phi_{1}/\delta+\mu)}}\mathbb{U}_{kx}$
$\displaystyle\mathbb{U}{}^{\dagger}K\mathbb{U}$ $\displaystyle=\epsilon_{k}.$
(52)

We emphasize that these results rely on many simplifications offered by a
constant spacetime offset $\phi=i\phi_{1}$. However, nearly identical
simplifications provide a similar evaluation for configurations with a
different constant field on each temporal slice.

Rather than compute derivatives of the action expressed with
$\log\det(\mathbbm{1}+\mathbb{F}_{\pm})$ (34) one may directly study $\log\det
M_{\pm}$ instead. In the case of constant field we may diagonalize $M$ with a
straightforward unitary Matsubara decomposition

$\displaystyle\mathbb{\Lambda}_{kn,xt}$
$\displaystyle=\mathbb{U}_{kx}e^{i\tilde{\omega}_{n}t}/\sqrt{N_{t}}$
$\displaystyle\tilde{\omega}_{n}$ $\displaystyle=(2n+1)\pi/N_{t}.$ (53)

One finds a structurally similar and numerically equal expression for the
Hessian in terms of matrix elements $T$

$\displaystyle\mathbb{H}_{x^{\prime}t^{\prime},xt}$
$\displaystyle=\left(\frac{1}{\tilde{U}}-1\right)\delta_{x^{\prime},x}\delta_{t^{\prime},t}-T_{+;xt,x^{\prime}t^{\prime}}T_{+;x^{\prime}t^{\prime},xt}-T_{-;xt,x^{\prime}t^{\prime}}T_{-;x^{\prime}t^{\prime},xt}$
(54) $\displaystyle T_{\pm;x^{\prime}t^{\prime},xt}$
$\displaystyle=\sum_{kn}\mathbb{\Lambda}{}^{\dagger}_{x^{\prime}t^{\prime},kn}\frac{e^{\pm(\delta\epsilon_{k}+\delta\mu+\phi_{1}+i\tilde{\omega}_{n})}}{1-e^{\pm(\delta\epsilon_{k}+\delta\mu+\phi_{1}+i\tilde{\omega}_{n})}}\mathbb{\Lambda}_{kn,xt}.$
(55)

We have numerically verified that these two formulations yield the same
Hessian.

## References

  * Splittorff and Verbaarschot [2007] K. Splittorff and J. Verbaarschot, Phase of the Fermion Determinant at Nonzero Chemical Potential, Phys. Rev. Lett. 98, 031601 (2007), arXiv:hep-lat/0609076 .
  * Schneider  _et al._ [2021] M. Schneider, J. Ostmeyer, K. Jansen, T. Luu, and C. Urbach, Simulating both parity sectors of the hubbard model with tensor networks, Physical Review B 104, 10.1103/physrevb.104.155118 (2021).
  * Kogut and Sinclair [2019] J. B. Kogut and D. K. Sinclair, Applying complex langevin simulations to lattice qcd at finite density, Phys. Rev. D 100, 054512 (2019), 1903.02622v3 .
  * Berger  _et al._ [2021] C. E. Berger, L. Rammelmüller, A. C. Loheac, F. Ehmann, J. Braun, and J. E. Drut, Complex langevin and other approaches to the sign problem in quantum many-body physics, Physics Reports 892, 1 (2021).
  * Cristoforetti  _et al._ [2013] M. Cristoforetti, F. Di Renzo, A. Mukherjee, and L. Scorzato, Monte Carlo simulations on the Lefschetz thimble: Taming the sign problem, Phys. Rev. D88, 051501 (2013), arXiv:1303.7204 [hep-lat] .
  * Cristoforetti  _et al._ [2014] M. Cristoforetti, F. Di Renzo, G. Eruzzi, A. Mukherjee, C. Schmidt, L. Scorzato, and C. Torrero, An efficient method to compute the residual phase on a Lefschetz thimble, Phys. Rev. D89, 114505 (2014), arXiv:1403.5637 [hep-lat] .
  * Lawrence [2018] S. Lawrence, Beyond Thimbles: Sign-Optimized Manifolds for Finite Density, PoS LATTICE2018, 149 (2018), arXiv:1810.06529 [hep-lat] .
  * Warrington [2019] N. C. Warrington,  _Taming the Sign Problem in Lattice Field Theory With Deformed Path Integral Contours_ , Ph.D. thesis, Maryland U. (2019).
  * Wynen  _et al._ [2019] J.-L. Wynen, E. Berkowitz, C. Körber, T. A. Lähde, and T. Luu, Avoiding ergodicity problems in lattice discretizations of the hubbard model, Physical Review B 100, 10.1103/physrevb.100.075141 (2019).
  * Detmold  _et al._ [2020] W. Detmold, G. Kanwar, M. L. Wagman, and N. C. Warrington, Path integral contour deformations for noisy observables, Phys. Rev. D 102, 014514 (2020), arXiv:2003.05914 [hep-lat] .
  * Alexandru  _et al._ [2020] A. Alexandru, G. Basar, P. F. Bedaque, and N. C. Warrington, Complex paths around the sign problem, arXiv preprint arXiv:2007.05436 (2020).
  * Detmold  _et al._ [2021] W. Detmold, G. Kanwar, H. Lamm, M. L. Wagman, and N. C. Warrington, Path integral contour deformations for observables in $SU(N)$ gauge theory, Phys. Rev. D 103, 094517 (2021), arXiv:2101.12668 [hep-lat] .
  * Rodekamp  _et al._ [2022] M. Rodekamp, E. Berkowitz, C. Gäntgen, S. Krieg, T. Luu, and J. Ostmeyer, Mitigating the hubbard sign problem with complex-valued neural networks, Phys. Rev. B 106, 125139 (2022), 2203.00390v2 .
  * Wynen  _et al._ [2021] J.-L. Wynen, E. Berkowitz, S. Krieg, T. Luu, and J. Ostmeyer, Leveraging machine learning to alleviate hubbard model sign problems, Phys. Rev. B 103, 125153 (2021), 2006.11221v1 .
  * Alexandru  _et al._ [2017] A. Alexandru, G. Basar, P. F. Bedaque, G. W. Ridgway, and N. C. Warrington, Monte carlo calculations of the finite density thirring model, Physical Review D 95, 10.1103/physrevd.95.014502 (2017), 1609.01730v2 .
  * Hubbard [1959] J. Hubbard, Calculation of partition functions, Physical Review Letters 3, 77 (1959).
  * Hubbard [1963] J. Hubbard, Electron correlations in narrow energy bands, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 276, 238 (1963).
  * Hubbard [1964a] J. Hubbard, Electron correlations in narrow energy bands. ii. the degenerate band case, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 277, 237 (1964a).
  * Hubbard [1964b] J. Hubbard, Electron correlations in narrow energy bands. iii. an improved solution, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 281, 401 (1964b).
  * Hubbard [1965] J. Hubbard, Electron correlations in narrow energy bands. iv. the atomic representation, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 285, 542 (1965).
  * Hubbard [1967a] J. Hubbard, Electron correlations in narrow energy bands. v. a perturbation expansion about the atomic limit, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 296, 82 (1967a).
  * Hubbard [1967b] J. Hubbard, Electron correlations in narrow energy bands. vi. the connexion with many-body perturbation theory, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 296, 100 (1967b).
  * Brower  _et al._ [2011] R. Brower, C. Rebbi, and D. Schaich, Hybrid Monte Carlo Simulation of Graphene on the Hexagonal Lattice (2011), arXiv:1101.5131 [hep-lat] .
  * Duane  _et al._ [1987] S. Duane, A. Kennedy, B. J. Pendleton, and D. Roweth, Hybrid monte carlo, Physics Letters B 195, 216 (1987).
  * Beyl  _et al._ [2018] S. Beyl, F. Goth, and F. F. Assaad, Revisiting the Hybrid Quantum Monte Carlo Method for Hubbard and Electron-Phonon Models, Phys. Rev. B 97, 085144 (2018), arXiv:1708.03661 [cond-mat.str-el] .
  * Ulybyshev  _et al._ [2020] M. Ulybyshev, C. Winterowd, and S. Zafeiropoulos, Lefschetz thimbles decomposition for the Hubbard model on the hexagonal lattice, Phys. Rev. D 101, 014508 (2020), arXiv:1906.07678 [cond-mat.str-el] .
  * Ulybyshev  _et al._ [2023] M. Ulybyshev, C. Winterowd, F. Assaad, and S. Zafeiropoulos, Instanton gas approach to the Hubbard model, Phys. Rev. B 107, 045143 (2023), arXiv:2207.06297 [cond-mat.str-el] .
  * Assaad and Herbut [2013] F. F. Assaad and I. F. Herbut, Pinning the order: the nature of quantum criticality in the Hubbard model on honeycomb lattice, Phys. Rev. X 3, 031010 (2013), arXiv:1304.6340 [cond-mat.str-el] .
  * Parisen Toldin  _et al._ [2015] F. Parisen Toldin, M. Hohenadler, F. F. Assaad, and I. F. Herbut, Fermionic quantum criticality in honeycomb and $\pi$-flux Hubbard models: Finite-size scaling of renormalization-group-invariant observables from quantum Monte Carlo, Phys. Rev. B 91, 165108 (2015), arXiv:1411.2502 [cond-mat.str-el] .
  * Otsuka  _et al._ [2016] Y. Otsuka, S. Yunoki, and S. Sorella, Universal quantum criticality in the metal-insulator transition of two-dimensional interacting Dirac electrons, Phys. Rev. X 6, 011029 (2016), arXiv:1510.08593 [cond-mat.str-el] .
  * Buividovich  _et al._ [2019] P. Buividovich, D. Smith, M. Ulybyshev, and L. von Smekal, Numerical evidence of conformal phase transition in graphene with long-range interactions, Phys. Rev. B 99, 205434 (2019), arXiv:1812.06435 [cond-mat.str-el] .
  * Buividovich  _et al._ [2018] P. Buividovich, D. Smith, M. Ulybyshev, and L. von Smekal, Hybrid-Monte-Carlo study of competing order in the extended fermionic Hubbard model on the hexagonal lattice, Phys. Rev. B 98, 235129 (2018), arXiv:1807.07025 [cond-mat.str-el] .
  * Ostmeyer  _et al._ [2020] J. Ostmeyer, E. Berkowitz, S. Krieg, T. A. Lähde, T. Luu, and C. Urbach, Semimetal–Mott insulator quantum phase transition of the Hubbard model on the honeycomb lattice, Phys. Rev. B 102, 245105 (2020), arXiv:2005.11112 [cond-mat.str-el] .
  * Ostmeyer  _et al._ [2021] J. Ostmeyer, E. Berkowitz, S. Krieg, T. A. Lähde, T. Luu, and C. Urbach, Antiferromagnetic character of the quantum phase transition in the Hubbard model on the honeycomb lattice, Phys. Rev. B 104, 155142 (2021), arXiv:2105.06936 [cond-mat.str-el] .
  * Saito  _et al._ [1998] R. Saito, G. Dresselhaus, and M. S. Dresselhaus,  _Physical Properties of Carbon Nanotubes_ (PUBLISHED BY IMPERIAL COLLEGE PRESS AND DISTRIBUTED BY WORLD SCIENTIFIC PUBLISHING CO., 1998) https://www.worldscientific.com/doi/pdf/10.1142/p080 .
  * Castro Neto  _et al._ [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, Rev. Mod. Phys. 81, 109 (2009), arXiv:0709.1163 [cond-mat.other] .
  * Alexandru  _et al._ [2018] A. Alexandru, P. F. Bedaque, H. Lamm, and S. Lawrence, Finite-Density Monte Carlo Calculations on Sign-Optimized Manifolds, Phys. Rev. D 97, 094510 (2018), arXiv:1804.00697 [hep-lat] .