# Bosonization of Majorana modes and edge states

Arkadiusz Bochniak<EMAIL_ADDRESS>Błażej Ruba
<EMAIL_ADDRESS>Jacek Wosiek<EMAIL_ADDRESS>Institute
of Theoretical Physics, Jagiellonian University, Poland.

###### Abstract

We present a bosonization procedure which replaces fermions with generalized
spin variables subject to local constraints. It requires that the number of
Majorana modes per lattice site matches the coordination number modulo two. If
this condition is not obeyed, then bosonization introduces additional
fermionic excitations not present in the original model. In the case of one
Majorana mode per site on a honeycomb lattice, we recover a sector of Kitaev’s
model. We discuss also decagonal and rectangular geometries and present
bosonization of the Hubbard model. For geometries with a boundary we find that
certain fermionic edge modes naturally emerge. They are of different nature
than edge modes encountered in topological phases of matter. Euclidean
representation for the unconstrained version of a spin system of the type
arising in our construction is derived and briefly studied by computing some
exact averages for small volumes.

††preprint: APS/123-QED

##  I Introduction

Bosonization is an old subject, which is of interests both in condensed matter
and high energy physics. The Jordan-Wigner transformation [1] is one of the
most famous methods. It provides an effective bosonization procedure in
$(1+1)$-dimensions. Its non-local character in higher dimensions leads to the
search of alternative methods. There exists a zoo of proposals [2, 3, 4, 5, 6,
7, 8, 9, 10, 11, 12, 13, 14, 15], including approaches motivated by the
Tomonaga-Luttinger model [16], generalizations or modifications of Witten’s
non-abelian bosonization [17], as well as purely algebraic approaches [18].
Bosonization is also closely related to the subject of dualities, such as the
Kramers-Wannier duality [19, 20] or the more recent web of dualities [11, 21,
22, 23].

Besides classical applications such as solving certain many body quantum
models exactly [24, 25] or overcoming sign problems in Monte Carlo studies
[26, 27, 28], bosonization has been invoked in the study of problems in
quantum computation [29, 30] and topological phases of matter [31, 32, 33, 34,
35].

In [4] a bosonization technique, here referred to as the $\Gamma$ model, was
proposed. It transforms in a local way a fermionic model into a generalized
spin system subject to constraints. This correspondence was then made more
precise in [36]. Generalization and new proofs were given in [5, 37].
Constraints present in the $\Gamma$ model were interpreted as the pure gauge
condition for a certain $\mathbb{Z}_{2}$ gauge field. Modification of these
constraints turned out to be equivalent to coupling fermions to an external
gauge field.

The most general version of the $\Gamma$ model developed so far is subject to
several important limitations. First, it corresponds to a system with one
fermion (hence two states) per lattice site. In this work we lift this
restriction and bosonize systems with arbitrary, not necessarily even, number
of Majorana modes (which are in a certain sense halves of an ordinary fermion)
per lattice site. Second, in the formulations considered before the present
work it was crucial that all lattice sites have an even number of neighbours.
This covers many interesting examples, inlcuding the toroidal geometries
frequently used in lattice simulations. Nevertheless, already for finite
square lattices (say, with open boundary conditions) there exists an issue
related to the existence of the boundaries. As remarked in [5, Appendix B],
coordination number changes for vertices on the boundary, which may call for
an adjustment of the bosonization procedure. It was argued that some Majorana
modes may be present on the boundary. In this paper we come back to this issue
and discuss it in detail.

We emphasize that the notion of a Majorana fermion used here has almost
nothing to do with the one from high energy physics [38], which refers to a
spinor field invariant under charge conjugation transformation. In particular
Lorentz symmetry (or lack thereof) plays no role. Here Majorana fermions are
self-adjoint operators obeying canonical anticommutation relations. Every
standard fermion may be decomposed into a pair of Majoranas (real and
imaginary part) in a canonical way; on the other hand pairing of Majoranas
into usual fermions depends on a choice of additional structure in the space
of Majorana modes 111More precisely, complex structure compatible with the
bilinear form determining the canonical anticommutation relations..

One source of interest in Majorana modes in physics comes from the theory of
superconductivity [40, 41]. In the presence of Abrikosov vortices [42] there
may exist a finite number of Majorana zero modes per vortex, which resemble
properties of Majorana particles [43]. Such phenomenon exists for eample for
chiral two-dimensional $p$-wave superconductors [44]. Analogous vortex-related
modes can be also found in superfluid ${}^{3}\mathrm{He}$ [45]. Majorana zero
modes are expected to appear also in the Moore–Read quantum Hall state [46]
with filling fraction $\nu=\frac{5}{2}$ (the so-called Pfaffian state). Quite
generally, free fermion systems characterized by nonzero $\mathbb{Z}_{2}$
topological invariant, such as the Kitaev’s quantum wires [47], are expected
to host Majorana zero modes on the boundary. Another exciting features of
Majorana modes is their potential in topological quantum computation [29, 48],
related to the possibility to realize non-abelian anyons. Feasibility of such
topological quantum computation is still being investigated [49].

The main idea underlying the $\Gamma$ model is to construct a representation
of the even subalgebra of fermionic operators (i.e. the subalgebra generated
by all bilinears) in terms of “spins” of sufficiently high dimension, or more
precisely in terms of Euclidean $\Gamma$ matrices (satisfying anticommutations
relations on-site, but otherwise commuting). As observed in [4], hopping
operators for fermions may be constructed given one $\Gamma(x,e)$ matrix per
lattice site $x$ for every edge $e$ incident to the given site. In addition,
one has to impose a certain constraints on states on the spin side. To
represent the standard algebra of fermions one has to specify, besides hopping
operators, also the fermionic parity operator on each site $x$. This operator
has to square to $1$ and anticommute with hopping operators along all edges
incident to $x$. In other words, one needs an additional $\Gamma$ matrix. If
$x$ has an even number of neighbours, this additional $\Gamma$ matrix may be
obtained (up to a trivial phase factor) simply by taking the product of all
$\Gamma(x,e)$ with fixed $x$. We emphasize that this construction does not
work if $x$ has an odd number of neighbours, since then the product of all
$\Gamma(x,e)$ commutes, rather than anticommutes with individual
$\Gamma(x,e)$. On the other hand introducing the additional $\Gamma$ matrix as
an independent object would lead to spurious degrees of freedom. Hence in this
version one restricts to even coordination numbers.

Generalization presented in this paper is based on a few simple observations.
First, in a system featuring an odd number of Majorana modes per lattice site
the on-site parity operators do not exist. Therefore the $\Gamma$ model on a
lattice with sites of odd degree should feature unpaired Majorana modes.
Second, in presence of multiple fermionic modes per site (say, due to spin or
orbital degeneracy) there exist additional independent bilinear operators
which can still be bosonized if one further increases the number of $\Gamma$
matrices per lattice site. In this way one obtains a mapping between fermions
and spins with only one requirement: the number of Majorana modes per site $x$
should be congruent modulo two to the number of neighbours of $x$. Even this
condition can be eventually lifted. Indeed, if it is not satisfied, one can
identify operators corresponding to spurious degrees of freedom and choose for
them trivial dynamics decoupled from the rest of the system.

It turns out that the $\Gamma$ model, in particular its version for arbitrary
lattices proposed in this paper, shares some features with models considered
in [50, 51] and [52], despite the fact that its origin and motivation were
different. We will now present a short comparison between these models. In
[50] the vector exchange model defined in terms of bond algebra was proposed.
Similarly like in our case and the models discussed by Kapustin et al. [2, 3],
the main idea was to say that two models (i.e. the fermionic and bosonic ones)
are equivalent if and only if their operator algebras are isomorphic. That is,
the statement was purely kinematic and Hamiltonian-independent. The idea of
using higher dimensional representations of Clifford algebras instead of the
Pauli matrices was introduced in [50] in order to define higher spin (e.g.
$\frac{3}{2}$, $\frac{7}{2}$ etc.) analogues of the Kiteav’s model. To proceed
with such fermionization procedure the need for lattices of coordination
number different than $3$ emerged. The interplay between the dimension of the
representation and the valency of lattice vertices is also discussed therein.
Relation between constraints and the choice of a $\mathbb{Z}_{2}$ gauge field
is also discussed and the counting of degrees of freedom is performed. In
contrast, we start with the fermionic theory and perform the bosonization
procedure based on the modification of the original $\Gamma$ model. The
(sector of) higher spin Kiteav’s model is a result of this procedure. We also
allow for multiple Majorana modes on different sites and this number may in
principle vary from site to site. As a consequence of the general bosonization
procedure, the relation of “unpaired” Majorana modes and generators of the
Clifford algebra associated to vertices is established. The fermionization
method analogous to the one in [50] was also proposed, at the same time, in
[51] and [52]. In the former case the periodic boundary conditions were
assumed, so that the role of the analogues of Polyakov lines discussed also in
details in [37] began to be important. The role of constraints was discussed,
together with the flux-attachment mechanism [53] and the interpretation of
modifying the constraints as a coupling to some external $\mathbb{Z}_{2}$
fields. We also remark that it was argued in [51] that models with a
$\mathbb{Z}_{2}$ gauge field chosen as in [5, Appendix B] may play a role for
$p$-wave superconductors. The discussion at the beginning of [52] is in the
same spirit as in [51]. In [52] the bulk-boundary correspondence is discussed
in more detail for such models. As pointed out in [50, 51, 52], these so-
called $\Gamma$-matrix models may have applications for spin liquids,
$(3+1)$-dimensional topological insulators and the $B$-phase of
${}^{3}\mathrm{He}$.

The organization of the paper is as follows. Details of our construction are
presented in Section II. Then in Section III we present examples: relation to
Kitaev’s model on hexagonal lattice, bosonization on a decagonal lattice and
bosonization of the Hubbard model on a rectangular lattice. Afterward, in
Section IV, we discuss boundary effects in the $\Gamma$ model. We describe the
example of square lattice with open boundary conditions and compare edge modes
identified there with those arising on the boundary of some topological phases
of matter.

In Section V an Euclidean representation of the simplest, unconstrained
$\Gamma$ model on a regular honeycomb lattice is proposed and briefly studied.
The time evolution generated by spin Hamiltonians considered here is more
complicated than in the standard Ising-like cases. Accordingly, Euclidean
three-dimensional spin systems emerging in this Section are unknown and
interesting by themselves. The feasibility of the standard, intermediate-
volume, Monte Carlo studies is crudely assessed on the basis of the exact
small-volume calculations.

##  II The bosonization method

We consider a lattice system with fermionic degrees of freedom, whose number
may vary from site to site. Real (Majorana) fermionic operators on the lattice
site $x$ will be denoted by $\psi_{\alpha}(x)$, with the index $\alpha$
running from $0$ to $n(x)$ with $n(x)\geq 0$. They are hermitian and satisfy
anticommutation relations

$\psi_{\alpha}(x)\psi_{\beta}(y)+\psi_{\beta}(y)\psi_{\alpha}(x)=2\delta_{x,y}\delta_{\alpha,\beta}.$
(1)

The total number of independent Majorana operators

$n=\sum_{x}\left(n(x)+1\right)$ (2)

is assumed to be even, which guarantees that

$(-1)^{F}=i^{\frac{n}{2}}\prod_{x}\prod_{\alpha=0}^{n(x)}\psi_{\alpha}(x)$ (3)

anticommutes with every fermion. We assume that $(-1)^{F}$ is an exactly
conserved quantity. Otherwise the Hamiltonian may be arbitrary.

We will be interested in the algebra of even operators (i.e. operators
commuting with $(-1)^{F}$). Any even operator may be expressed as a linear
combination of products of bilinears of the following two types:

$\displaystyle S(e)$ $\displaystyle=\psi_{0}(x)\psi_{0}(y)\text{ for an edge
}e\text{ from }x\text{ to }y,$ (4a) $\displaystyle T_{\alpha}(x)$
$\displaystyle=\psi_{0}(x)\psi_{\alpha}(x)\text{ for }\alpha\neq 0.$ (4b)

All $S$ and $T$ operators are skew-hermitian and square to $-1$. Furthermore
we have that

  * •

$S(e)S(e^{\prime})=\pm S(e^{\prime})S(e)$, with the minus sign only if $e$
shares exactly one endpoint with $e^{\prime}$,

  * •

$S(e)T_{\alpha}(x)=\pm T_{\alpha}(x)S(e)$, with the minus sign only if $x$ is
incident to $e$,

  * •

$T_{\alpha}(x)T_{\beta}(y)=\pm T_{\beta}(y)T_{\alpha}(x)$, with the minus sign
only if $x=y$ and $\alpha\neq\beta$.

It can be shown [5] that all relations in the algebra of even operators are
generated by those given above and what will be called loop relations: if
edges $e_{1},\ldots,e_{m}$ form a loop (i.e. $e_{i}$ terminates at the initial
point of $e_{i+1}$, with the convention that $e_{m+1}=e_{1}$), then

$S(e_{1})\ldots S(e_{m})=1.$ (5)

To bosonize the system, we generalize the approach proposed in [5]. For each
lattice site $x$ we construct a Clifford algebra with generators
$\Gamma(x,e)$, one for each edge $e$ incident to $x$, and
$\Gamma^{\prime}_{\alpha}(x)$ with $0<\alpha\leq n(x)$. They are hermitian
matrices satisfying

$\displaystyle\Gamma(x,e)\Gamma(x,e^{\prime})+\Gamma(x,e^{\prime})\Gamma(x,e)$
$\displaystyle=2\delta_{e,e^{\prime}},$ (6a)
$\displaystyle\Gamma^{\prime}_{\alpha}(x)\Gamma^{\prime}_{\beta}(x)+\Gamma^{\prime}_{\beta}(x)\Gamma^{\prime}_{\alpha}(x)$
$\displaystyle=2\delta_{\alpha,\beta},$ (6b)
$\displaystyle\Gamma(x,e)\Gamma^{\prime}_{\alpha}(x)+\Gamma^{\prime}_{\alpha}(x)\Gamma(x,e)$
$\displaystyle=0.$ (6c)

Gamma matrices located on distinct lattice sites are taken to commute, and the
full Hilbert space is the tensor product of on-site Hilbert spaces. In this
sense the new system is bosonic. Fermionic bilinears are mapped to bosonic
operators according to the local prescription

$\displaystyle\widehat{S}(e)=i\,\Gamma(x,e)\Gamma(y,e)\text{ for an edge
}e\text{ from }x\text{ to }y,$ (7a)
$\displaystyle\widehat{T}_{\alpha}(x)=i\,\Gamma^{\prime}_{\alpha}(x)\text{ for
}\alpha\neq 0,$ (7b)

where the hat serves as an indicator that we are referring to the bosonized
operators, rather than those in the original fermionic system. It is
straightforward to check that $\widehat{S}$ and $\widehat{T}$ operators
satisfy all relations obeyed by $S$ and $T$, except for the loop relations.
Instead, for a loop $\ell$ formed by edges $e_{1},\ldots,e_{m}$ one has that
the operator

$W(\ell)=\widehat{S}(e_{1})\ldots\widehat{S}(e_{m})$ (8)

squares to $1$ and commutes with all $\widehat{S}$ and $\widehat{T}$. We are
forced to impose the constraint

$W(\ell)|\mathrm{phys}\rangle=|\mathrm{phys}\rangle\text{ for every loop
}\ell.$ (9)

We remark that modifying the constraint to the form

$W(\ell)|\mathrm{phys}\rangle=\omega(\ell)|\mathrm{phys}\rangle$ (10)

with prescribed $\omega(\ell)=\pm 1$ is equivalent [5] to coupling fermions to
a background $\mathbb{Z}_{2}$ gauge field for the $(-1)^{F}$ symmetry, such
that $\omega(\ell)$ is the holonomy along $\ell$.

Now let $\deg(x)$ be the number of neighbors of a site $x$ and put
$N(x)=\deg(x)+n(x)$. We consider the operator

$\gamma(x)=i^{\frac{N(x)(N(x)-1)}{2}}\prod_{e}\Gamma(x,e)\prod_{\alpha\neq
0}\Gamma^{\prime}_{\alpha}(x).$ (11)

Its phase is chosen so that $\gamma(x)^{2}=1$. If $N(x)$ is odd, $\gamma(x)$
commutes with all gamma matrices, so one may impose a relation $\gamma(x)=1$
or $\gamma(x)=-1$. This amounts to choosing one of two irreducible
representation of the Clifford algebra on $x$. If $N(x)$ is even, $\gamma(x)$
anticommutes with all gamma matrices, so it defines an additional gamma
matrix: $\gamma(x)=\Gamma^{\prime}_{n(x)+1}(x)$. In this case Eqs. (7) provide
a bosonization of a system featuring one more Majorana fermion on the site $x$
than we have had originally. Therefore formally we bosonize only systems with
all $N(x)$ odd, but the case in which this condition is not satisfied may be
handled by choosing for the spurious fermions a trivially gapped Hamiltonian,
not interacting with the original fermions.

Bosonic system with constraints imposed is equivalent to the sector of the
fermionic system (possibly including the spurious fermions) characterized by
one of the two possible values of $(-1)^{F}$, defined including the spurious
fermions. Which possibility is realized depends on the lattice geometry and
the choice of values of $\gamma$ operators. In the remainder of this section
we sketch the proof of this fact, while details have been given in [5] in a
slightly less general context.

First, on the space of solutions of the constraints all relations satisfied by
$S$ and $T$ operators are obeyed by $\widehat{S}$ and $\widehat{T}$. Therefore
this space is a representation of the algebra of even operators. Every
representation of this algebra is a direct sum of irreducible representations,
which are the two halves of the Fock space described by two values of
$(-1)^{F}$. We will argue that only one of the two irreducible representations
actually occurs in the decomposition and that the multiplicity is equal to
one.

For the first part of the claim, it is sufficient to observe that the product
of all $S$ and $T$ operators is proportional to the fermionic parity operator,
while the product of all $\widehat{S}$ and $\widehat{T}$ is proportional to
the product of all gamma matrices. The latter is proportional to $1$, because
for every lattice site we have $\gamma(x)=1$ or $\gamma(x)=-1$. Combining
these two results we conclude that $(-1)^{F}$ is represented by a scalar
operator in the bosonic model. It is possible to determine whether it is equal
to $+1$ or $-1$ by tracking phases carefully in the above argument. Details
depend on the lattice geometry.

For the second part of the claim it suffices to calculate the dimension of the
space of solutions of constraints. This is facilitated by considering also
modified forms of constraints. Those are in one-to-one correspondence with
gauge orbits of background $\mathbb{Z}_{2}$ gauge fields. Thus there are
$2^{N_{1}-N_{0}+1}$ of them, where $N_{0}$ is the number of lattice sites and
$N_{1}$ is the number of edges. One can show that spaces of solutions of
modified constraints all have the same dimension by explicitly constructing
unitary operators which map between them. Therefore the dimension of the space
of solutions of constraints (9) is equal to the dimension of the whole Hilbert
space in the bosonic model divided by $2^{N_{1}-N_{0}+1}$. Using the well-
known values of dimensions of irreducible representations of Clifford algebras
we find that there are $2^{\frac{n}{2}-1}$ linearly independent solutions of
constraints. This number is equal to the dimension of one half of the Fock
space.

##  III Examples

We will now show how the general construction presented in the previous
section works in specific examples. We begin with a model defined on the
honeycomb lattice and discuss its relation with the Kitaev’s model [54]. Then
we discuss its three-dimensional deformation, the decagonal lattice.
Discussion of boundary effects is postponed to Section IV.

We stress that this choice of lattices has been made mostly in order to
simplify the presentation. Bosonization prescription from Sec. II is valid
also on lattices of more complicated geometry: coordination numbers may vary
from site to site and the translation symmetry is not necessary.

###  III.1 Honeycomb lattice and Kitaev’s model

In this subsection we present an example involving a honeycomb lattice. It is
arguably the simplest two dimensional lattice with vertices of odd degree.
Since our construction is sensitive only to the topology rather than geometry
of the lattice, honeycomb lattice is equivalent to the brick wall lattice, see
Fig. 1.

(a)

(b)

Figure 1: (a) Honeycomb and (b) brick wall lattices are geometrically
different, but topologically equivalent.

We will bosonize a system featuring one Majorana fermion $\psi$ per lattice
site. This requires three gamma matrices. They can be represented by Pauli
matrices $\sigma_{X},\sigma_{Y},\sigma_{Z}$, which are assigned to edges of
the lattice as illustrated in Fig. 2. Thus for a lattice site $x$ and an edge
$e$ labeled by $I\in\\{X,Y,Z\\}$ we have $\Gamma(x,e)=\sigma_{I}(x)$.

Figure 2: The assignment of Pauli matrices.

For a plaquette $P$ depicted in Fig. 3, the corresponding constraint takes the
form $W_{P}|\mathrm{phys}\rangle=-|\mathrm{phys}\rangle$, where

$W_{P}=\sigma_{X}(x_{1})\sigma_{Y}(x_{2})\sigma_{Z}(x_{3})\sigma_{X}(x_{4})\sigma_{Y}(x_{5})\sigma_{Z}(x_{6}).$
(12)

Those are the Kitaev’s plaquette operators [54].

Figure 3: A plaquette $P$ of the honeycomb lattice.

As a specific example, let us consider the Hamiltonian

$H=i\sum\limits_{I\in\\{X,Y,Z\\}}\sum_{\begin{subarray}{c}\text{type }I\\\
\text{edges}\end{subarray}}J_{I}\,\psi(x)\psi(y),$ (13)

where $x,y$ are the endpoints of the given edge, and $J_{X},J_{Y}$ and $J_{Z}$
are parameters of the model. According to the prescription given in Eq. (7),
it corresponds to the spin Hamiltonian

$\widehat{H}=-\sum\limits_{I\in\\{X,Y,Z\\}}\sum_{\begin{subarray}{c}\text{type
}I\\\ \text{edges}\end{subarray}}J_{I}\,\sigma_{I}(x)\sigma_{I}(y),$ (14)

subject to the constraint $W_{P}=-1$ for every plaquette $P$. This Hamiltonian
has been proposed in [54], where its study was reduced to diagonalization of
quadratic fermionic Hamiltonians. Our approach provides an alternative
derivation of this result. Subspaces defined by different values of $W_{P}$
correspond to the Hamiltonian $H$ modified by including a background
$\mathbb{Z}_{2}$ gauge field.

###  III.2 Decagonal lattice

An example of a three-dimensional trivalent lattice is provided by the
decagonal geometry. A convenient representation is shown in Fig. 4, where one
layer of such lattice is presented, together with edges connecting it with the
adjacent layers. It can be thought of as a deformation of the brick wall
lattice. In the brick wall geometry, each red site was connected with a green
one to its north, while in the decagonal geometry it is instead connected with
a green site lying in the layer underneath.

Figure 4: One layer of the decagonal lattice.

By the similarity with the brick wall geometry, one can easily generalize the
results from Section III.1. Using the identification of edges between
honeycomb and brick wall lattices we attach Pauli matrices to pairs $(x,e)$ of
the decagonal lattice, see Fig. 5.

(a)

(b)

Figure 5: The assignment of Pauli matrices for (a) the brick wall lattice and
(b) the decagonal lattice.

As an example, constraint associated to the plaquette from Fig. 6 takes the
form $W_{P}|\mathrm{phys}\rangle=-|\mathrm{phys}\rangle$, where

$\displaystyle W_{P}=$
$\displaystyle\sigma_{X}(x_{1})\sigma_{Y}(x_{2})\sigma_{Z}(x_{3})\sigma_{Z}(x_{4})\sigma_{Z}(x_{5})$
(15)
$\displaystyle\times\sigma_{X}(x_{6})\sigma_{Y}(x_{7})\sigma_{Z}(x_{8})\sigma_{Z}(x_{9})\sigma_{Z}(x_{10}).$

There exist also plaquettes not contained within one layer, but for the sake
of brevity we will not write down the explicit formulas.

Figure 6: Plaquette $P$ within a single layer of the decagonal lattice.

As in the honeycomb lattice case, every edge is labeled by $I\in\\{X,Y,Z\\}$
and there is a correspondence between the Hamiltonians in Eqs. (13) and (14).

###  III.3 Hubbard model

In the preceding examples only one kind of fermionic variables was involved.
Here we discuss the simplest model with an additional quantum number involved
- the Hubbard model [55] on the square lattice.

The Hamiltonian of the Hubbard model consists of two terms, $H=H_{0}+V$, where

$\displaystyle H_{0}=-t\sum\limits_{\langle
xy\rangle}\sum\limits_{\sigma=\uparrow,\downarrow}\left(c_{\sigma}^{\dagger}(x)c_{\sigma}(y)+c^{\dagger}_{\sigma}(y)c_{\sigma}(x)\right),$
(16a) $\displaystyle V=U\sum\limits_{x}n_{\uparrow}(x)n_{\downarrow}(x),$
(16b)

where the edge connecting sites $x$ and $y$ is denoted by $\langle xy\rangle$.
Here $c_{\sigma}^{\dagger}(x)$ creates a fermion with spin $\sigma$ at
position $x$, $n_{\sigma}(x)=c_{\sigma}^{\dagger}(x)c_{\sigma}(x)$, and
$t,U\in\mathbb{R}$. These fermionic operators may be decomposed into Majoranas
as

$\displaystyle c_{\uparrow}(x)$
$\displaystyle=2^{-1}(\psi_{0}(x)+i\psi_{1}(x)),$ (17a) $\displaystyle
c^{\dagger}_{\uparrow}(x)$ $\displaystyle=2^{-1}(\psi_{0}(x)-i\psi_{1}(x)),$
(17b) $\displaystyle c_{\downarrow}(x)$
$\displaystyle=2^{-1}(\psi_{2}(x)+i\psi_{3}(x)),$ (17c) $\displaystyle
c^{\dagger}_{\downarrow}(x)$ $\displaystyle=2^{-1}(\psi_{2}(x)-i\psi_{3}(x)).$
(17d)

To perform bosonization we need seven $\Gamma$ matrices per site, six of which
are independent: the seventh may be taken to be the product of the first six
and the imaginary unit. The bosonized Hamiltonian takes the form:

$\displaystyle\widehat{H}_{0}=-\frac{it}{2}\sum\limits_{e=\langle
xy\rangle}\Gamma(x,e)\Gamma(y,e)\left(\Gamma_{1}^{\prime}(x)-\Gamma_{1}^{\prime}(y)-i\Gamma_{2}^{\prime}(x)\Gamma_{3}^{\prime}(y)+i\Gamma_{3}^{\prime}(x)\Gamma_{2}^{\prime}(y)\right),$
(18a)
$\displaystyle\widehat{V}=\frac{U}{4}\sum\limits_{x}\left(1-\Gamma_{1}^{\prime}(x)\right)\left(1+i\Gamma_{2}^{\prime}(x)\Gamma_{3}^{\prime}(x)\right).$
(18b)

Figure 7: The assignment of gamma matrices. Figure 8: A plaquette $p$ of the
rectangular lattice.

To write down the constraints, it is covenient to denote $\Gamma(x,e)$ as
$\Gamma_{i}(x)$ for edge $i$ pointing from $x$ in the $i$-th direction,
$i\in\\{\pm 1,\pm 2\\}$ (see Figure 7). Now consider a plaquette $P$ as in
Fig. 8. The corresponding constraint takes the form
$W_{P}|\mathrm{phys}\rangle=-|\mathrm{phys}\rangle$, where

$W_{P}=\Gamma_{1,2}(x_{1})\Gamma_{2,-1}(x_{2})\Gamma_{-1,-2}(x_{3})\Gamma_{-2,1}(x_{4}).$
(19)

Here we abbreviated $\Gamma_{i,j}(x):=\Gamma_{i}(x)\Gamma_{j}(x)$. We note
that this constraint does not at all involve primed gamma matrices, which we
had to introduce in order to implement multiple fermions per site. It is
characteristic for the square lattice geometry.

One annoying feature of the presented construction is that spin up and spin
down states are not treated completely symmetrically. Nevertheless, symmetries
of the Hubbard model are implemented also in the bosonized model. First, we
have conservation of the total particle number. Particle number on a single
lattice site $x$ is bosonized in the following way:

$\sum_{\sigma}c_{\sigma}^{\dagger}c_{\sigma}\longleftrightarrow
1-\frac{1}{2}\Gamma_{1}^{\prime}+\frac{i}{2}\Gamma_{2}^{\prime}\Gamma_{3}^{\prime}.$
(20)

For spin ($\mathrm{SU}(2)$ generators) operators we have:

$\displaystyle
c_{\uparrow}^{\dagger}c_{\uparrow}-c_{\downarrow}^{\dagger}c_{\downarrow}\longleftrightarrow-\frac{1}{2}\Gamma_{1}^{\prime}-\frac{i}{2}\Gamma_{2}^{\prime}\Gamma_{3}^{\prime},$
(21a) $\displaystyle
c_{\uparrow}^{\dagger}c_{\downarrow}\longleftrightarrow\frac{i}{4}(1-\Gamma_{1}^{\prime})(\Gamma_{2}^{\prime}+i\Gamma_{3}^{\prime}),$
(21b) $\displaystyle
c_{\downarrow}^{\dagger}c_{\uparrow}\longleftrightarrow-\frac{i}{4}(1+\Gamma_{1}^{\prime})(\Gamma_{2}^{\prime}-i\Gamma_{3}^{\prime}).$
(21c)

It noteworthy that all on-site symmetries of the original fermionic model are
also on-site after bosonization. Furthermore, the corresponding charges are
expressed entirely in terms of _primed_ gamma matrices. This is a general
property of our construction.

##  IV Boundary effects

###  IV.1 Rectangular lattice with a boundary

We will now discuss bosonization of a system on an $L_{x}\times L_{y}$
rectangular lattice, with two Majorana fermions $\psi_{0}$, $\psi_{1}$ per
lattice site. Every site $x$ in the bulk has four neighbors, corresponding to
four gamma matrices $\Gamma_{\pm i}(x)$, $i=1,2$, as in our discussion of the
Hubbard model. According to the prescription given in Sec. II, we need also an
additional gamma matrix $\Gamma^{\prime}_{1}(x)$. It can be eliminated by
imposing relations discussed below equation (11). We choose the convention
$\Gamma_{1}^{\prime}(x)=\Gamma_{-1}(x)\Gamma_{1}(x)\Gamma_{-2}(x)\Gamma_{2}(x)$.
Therefore in the end we need only unprimed gamma matrices. Constraints are
identical as in the discussion of the Hubbard model.

More explicitly, our bosonization prescription reads

  * •

$i\Gamma_{1}(x)\Gamma_{-1}(y)\longleftrightarrow\psi_{0}(x)\psi_{0}(y)$ if $y$
is the eastern neighbor of $x$,

  * •

$i\Gamma_{2}(x)\Gamma_{-2}(y)\longleftrightarrow\psi_{0}(x)\psi_{0}(y)$ if $y$
is the northern neighbor of $x$,

  * •

$i\Gamma_{-1}(x)\Gamma_{1}(x)\Gamma_{-2}(x)\Gamma_{2}(x)\longleftrightarrow\psi_{0}(x)\psi_{1}(x)$.

Now, we look closely at the situation on the boundary. First, sites on the
southern edge (see Fig. 9) have no neighbors in the direction $-2$. We may
reinterpret the $\Gamma_{-2}$ matrix as an additional $\Gamma^{\prime}$,
corresponding to a spurious Majorana fermion on the boundary. More precisely,
for every site $x_{i}$ on the southern edge we introduce an additional
Majorana operator $\chi_{\mathrm{S}}(x_{i})$. Bosonization prescription for
$\chi_{S}$ fermions takes the form

$i\Gamma_{-2}(x_{i})\longleftrightarrow\psi_{0}(x_{i})\chi_{\mathrm{S}}(x_{i}).$
(22)

Similarly for the northern, eastern and western edges we introduce Majorana
fermions $\chi_{\mathrm{N}}$, $\chi_{\mathrm{E}}$ and $\chi_{\mathrm{W}}$.

Figure 9: The southern boundary of the rectangular lattice.

At each of the four corners (which are geometrically of codimension two) there
are two $\chi$ fermions. For example the south-east corner $x_{\mathrm{SE}}$
hosts four Majorana operators $\psi_{0}(x_{\mathrm{SE}})$,
$\psi_{1}(x_{\mathrm{SE}})$, $\chi_{S}(x_{\mathrm{SE}})$ and
$\chi_{E}(x_{\mathrm{SE}})$.

We now determine the identity resulting from existence of the boundary. First,
notice that for every lattice site
$x=(a,b)\in\\{1,\ldots,L_{x}\\}\times\\{1,\ldots,L_{y}\\}$ we have

$\prod\limits_{a=1}^{L_{x}-1}\psi_{0}(a,b)\psi_{0}(a+1,b)=\psi_{0}(a,b)\psi_{0}(L_{x},b).$
(23)

Consequently, our bosonization prescription yields

$\begin{split}&\psi_{0}(1,b)\psi_{0}(L_{x},b)\\\ &\longleftrightarrow
i^{L_{x}-1}\Gamma_{-1}(1,b)\left(\prod\limits_{a=1}^{L_{x}}\Gamma_{-1,1}(a,b)\right)\Gamma_{1}(L_{x},b)\end{split}$
(24)

and

$\begin{split}&\psi_{0}(a,1)\psi_{0}(a,L_{y})\\\ &\longleftrightarrow
i^{L_{y}-1}\Gamma_{-2}(a,1)\left(\prod\limits_{b=1}^{L_{y}}\Gamma_{-2,2}(a,b)\right)\Gamma_{1}(a,L_{y}),\end{split}$
(25)

for every $1\leq b\leq L_{y}$ and $1\leq a\leq L_{x}$, respectively. By a
straightforward computation one can check that it results in the following
correspondence

$\begin{split}&i^{(L_{x}-L_{y})^{2}+2(L_{x}+L_{y})}\prod\limits_{a,b}\Gamma_{-2,2,-1,1}(a,b)\\\
&\longleftrightarrow\chi_{\partial_{S}}\chi_{\partial_{N}}\chi_{\partial_{W}}\chi_{\partial_{E}},\end{split}$
(26)

where we have introduced the abbreviated notation
$\chi_{\partial_{S}}=\prod\limits_{a=1}^{L_{x}}\chi_{S}(a,1)$, and so on.
Since $\Gamma_{-2,2,-1,1}(x)$ corresponds to $-i\psi_{0}(x)\psi_{1}(x)$, which
is the parity operator $(-1)^{F_{\psi}(x)}$ for $\psi$ fermions at site $x$,
we end up with the following constraint

$(-1)^{F_{\psi}}=\kappa\chi_{\partial_{S}}\chi_{\partial_{N}}\chi_{\partial_{W}}\chi_{\partial_{E}},$
(27)

where $\kappa=i^{-(L_{x}-L_{y})^{2}+2(L_{x}+L_{y})}$ is a geometrical phase
factor. In particular, for lattices with $L_{x}\equiv L_{y}\pmod{2}$ we have
$\kappa=1$.

Summarizing, we started from the system of $\psi$ fermions, but our
bosonization gave us a bosonic system equivalent to $\psi$ fermions together
with $\chi$ fermions on the boundary. Since operators corresponding on the
spin side to spurious $\chi$ modes have been identified, this is not a
problem. Indeed, bosonizing suitable Hamiltonian for $\chi$ fermions using
formula (22) and its analogues for other components of the boundary we may
make $\chi$ arbitrarily heavy, e.g. dimerized with large dissociation energy.
This provides a physical interpretation for additional constraints on the
boundary introduced in [5, Appendix B]. After imposing them, one obtains a
spin system corresponding to $\psi$ fermions on a lattice with boundary.

###  IV.2 Comment about topological phases

One of the remarkable features of many topologically nontrivial phases of
matter is their interesting (robust) behaviour on the boundaries, which are of
dimension $d-1$. This is the case in particular for free fermion systems, for
which one has a well-established bulk-boundary correspondence: topological
invariants in the bulk signal existence of modes localized near the boundary,
responsible for closing the gap in finite volume.

Robustness of the boundary modes in often interpreted as manifestation of an
anomaly of the boundary theory. Presence of the anomaly implies existence of
some degrees of freedom “saturating” the anomaly. On the other hand, it is
also expected that the anomalous $(d-1)$-dimensional system is inconsistent on
its own: it may be realized only on the boundary of a $d$-dimensional system.
As an example, chiral fermions cannot be realized on the lattice (this is the
Nielsen-Ninomiya theorem, see [56, 57]), but they may exist on the boundary or
domain wall (more generally, a defect) in a higher dimensional system. This is
at the heart of the bulk boundary correspondence.

On the other hand, boundary modes found in our bosonization prescription
correspond to a standalone (hence “non-anomalous”) system on the boundary. Now
suppose that we bosonize a free fermion system with a nonzero topological
invariant, say on a half-space. Then on the boundary we will have boundary
modes predicted by the bulk-boundary correspondence and $\chi$ fermions
described in the previous subsection. We can gap out the latter fermions by
including in the Hamiltonian a suitable term localized in the boundary. This
is believed not to be true for the former. In the case of invariants which
remain robust in presence of interactions [58] it is natural to expect that
even after including a coupling between $\chi$ fermions and $\psi$ fermionis,
boundary modes originating from a topological invariant will persist.

Of course the full picture of topological invariants and boundary modes has to
involve the choice of a Hamiltonian, or at least some class of Hamiltonians.
On the other hand, the discussion presented here is mostly concerned with
properties of algebras of observables. It would be interesting to understand
better the relation between bosonization and bulk-boundary correspondence.
Such questions are relevant, for example, for the problem of discretization of
chiral fermions.

##  V Euclidean representation of unconstrained “Majorana spins”

The next goal is to construct an Euclidean Ising-like action, with two
different couplings, $\beta_{t}$ and $\beta_{s}$, which in the continuous time
limit

$\beta_{t}\rightarrow\infty,\quad\epsilon=e^{-\beta_{t}}\rightarrow
0,\quad\beta_{s}=\epsilon\lambda\rightarrow 0$ (28)

is described by the Hamiltonian (14) with parameters $J_{X}=J_{Y}=1$,
$J_{Z}=\lambda$ (we will not impose constraints at this point). Following [59,
60], this is done by demanding that the transfer matrix elements defined by
the Boltzmann weight

$\langle s^{\prime}|T|s\rangle=e^{-\mathsf{L}(s^{\prime},s)}$ (29)

coincide with these of the Euclidean evolution operator

$\langle s^{\prime}|T|s\rangle=\langle s^{\prime}|e^{-\epsilon
H}|s\rangle=\langle s^{\prime}|1-\epsilon H|s\rangle$ (30)

up to terms of order $\epsilon$. Here $\epsilon$ is the elementary time step
and $s$ and $s^{\prime}$ denote configurations of spins at subsequent time
instants.

The time evolution generated by (14) consists of elementary double-spin flips,
in contrast to the Ising system in which the dynamics is driven by single spin
flips. In order to gain some orientation in this problem, we start by deriving
an Euclidean action for a simpler, one-dimensional quantum Hamiltonian

$\begin{split}H_{1d}=&-\sum_{k}\sigma_{X}(x_{k})\sigma_{X}(x_{k+1})\\\
&-\lambda\sum_{k}\sigma_{Z}(x_{k})\sigma_{Z}(x_{k+1}).\end{split}$ (31)

###  V.1 Basic idea and the $(1+1)$-dimensional example

In the Ising model, the basic trick relating Hamiltonian and functional
formulations is to classify all variations of a multiple-spin state into
classes with fixed number of single spin flips.

On the Euclidean side, the number of single flips between two time slices is
counted by the two-row action

$\mathsf{L}_{1}(s^{\prime},s)=\frac{1}{2}\sum_{k}(1-s_{k}s_{k}^{\prime}).$
(32)

It corresponds to the Hamiltonian

$-\sum_{k}\sigma_{X}(x_{k}).$ (33)

In the present case (31), we are seeking to single out the double spin flips
out of all possible changes of a row of spins. Therefore we begin with the
Euclidean eight-spin action which counts isolated double flips

$\begin{split}\mathsf{L}_{2}^{(8)}&=\frac{1}{2^{4}}\sum_{k}(1+s_{k-1}s_{k-1}^{\prime})(1-s_{k}s_{k}^{\prime})\\\
&\times(1-s_{k+1}s_{k+1}^{\prime})(1+s_{k+2}s_{k+2}^{\prime}),\end{split}$
(34)

Simpler functions can be also used, hence we shall omit the “$(8)$”
superscript if not necessary.

We need to arrange the final, Euclidean action such that in the continuous
time limit it gives weight $\epsilon$ to double flips while all other, single
and multiple, flips are of higher order in $\epsilon=e^{-\beta_{t}}$. This is
achieved by the combination

$\mathsf{L}_{\mathrm{kin}}(s^{\prime},s)=\beta_{t}\left(p(\mathsf{L}_{1}-2\mathsf{L}_{2})+\mathsf{L}_{2}\right),$
(35)

where $p\geq 2$ is a free parameter.

It is easy to check that $\mathsf{L}_{2}^{(8)}$ may also be replaced in (35)
by the simpler function

$\mathsf{L}_{2}^{(6)}=\frac{1}{8}\sum_{k}(1+s_{k-1}s_{k-1}^{\prime})(1-s_{k}s_{k}^{\prime})(1-s_{k+1}s_{k+1}^{\prime}).$
(36)

This definition prescribes different weights to non-leading transitions, but
results in the same double flip kinetic part of (31). This is an illustration
of the well known fact that many different Euclidean discretizations have the
same continuous time limit, hence also the same Hamiltonian.

Action for a single transition has to be supplemented by a potential term:

$\displaystyle\mathsf{L}(s^{\prime},s)=\mathsf{L}_{\mathrm{kin}}(s^{\prime},s)+\mathsf{L}_{\mathrm{pot}}(s^{\prime},s),$
(37a)
$\displaystyle\mathsf{L}_{\mathrm{pot}}(s^{\prime},s)=-\frac{\beta_{s}}{2}\sum_{k}\left(s_{k}s_{k+1}+s_{k}^{\prime}s_{k+1}^{\prime}\right).$
(37b)

Complete Euclidean action for the $L_{x}\times L_{t}$ spins is obtained by
composing elementary transfer matrices, which amounts to adding the
corresponding actions:

$S(s(L_{t}),\ldots,s(1))=\sum\limits_{t=1}^{L_{t}}\mathsf{L}(s(t+1),s(t)).$
(38)

This concludes our construction of the two-dimensional, Euclidean system which
in the continuum time limit is described by the Hamiltonian (31).

####  V.1.1 $\sigma_{Y}\sigma_{Y}$ terms - the phases.

The second example deals with the phase generating kinetic terms

$\displaystyle H_{1d}^{\mathrm{ph}}=$
$\displaystyle-\sum_{k\;\mathrm{even}}\sigma_{X}(x_{k})\sigma_{X}(x_{k+1})-\sum_{k\;\mathrm{odd}}\sigma_{Y}(x_{k})\sigma_{Y}(x_{k+1})$
$\displaystyle-\lambda\sum_{k}\sigma_{Z}(x_{k})\sigma_{Z}(x_{k+1})$ (39)

still in one space dimension.

Begin with an evolution of a two spin system:

$s=\\{s_{1},s_{2}\\}\rightarrow
s^{\prime}=\\{s_{1}^{\prime},s_{2}^{\prime}\\}.$ (40)

As far as the change of spin states is considered, the action of
$\sigma_{Y}\sigma_{Y}$ is the same as that of $\sigma_{X}\sigma_{X}$. The only
difference is a phase factor:

$\begin{split}&\sigma_{Y}(x_{1})\sigma_{Y}(x_{2})|s_{1},s_{2}\rangle\\\
&=\exp{\left(\frac{i\pi}{2}(s_{1}+s_{2})\right)}\sigma_{X}(x_{1})\sigma_{X}(x_{2})|s_{1},s_{2}\rangle.\end{split}$
(41)

Generalization to a whole row of $L$ spins is straightforward. The kinetic
term of the Hamiltonian (39) will be reproduced by the action (35)
supplemented by a phase (41) for each odd edge. This gives for the new action
of the two complete rows (and with the unchanged diagonal potential term)

$\begin{split}\mathsf{L}_{\mathrm{ph}}(s^{\prime},s)&=\beta_{t}\left(2(S_{1}-2S_{2})+S_{2}\right)+\beta_{s}S^{\mathrm{pot}}\\\
&+\frac{i\pi}{2}\frac{1}{2^{3}}\sum_{x-\mathrm{odd}}(s_{x}+s_{x+1})(1+s_{x-1}s_{x-1}^{\prime})(1-s_{x}s_{x}^{\prime})(1-s_{x+1}s_{x+1}^{\prime}).\end{split}$
(42)

###  V.2 $(2+1)$-dimensional system

As a $(2+1)$-dimensional example we consider here the honeycomb lattice
discussed in III.1. As remarked therein, it is convenient to represent it in a
brick wall form, cf. Fig. 1.

Our Hamiltionian is of the form (14) with parameters $J_{X}=J_{Y}=1$ and
$J_{Z}=\lambda$, and can be written as a sum of two terms

$H=H_{\mathrm{kin}}+\lambda H_{\mathrm{pot}},$ (43)

where the kinetic term contains sums over all edges of type $X$ and $Y$, while
the potential one is a sum over edges of type $Z$. The labelling of the edges
of the brick wall lattice is shown explicitly in Fig. 10 and is consistent
with the one in Fig. 2.

Figure 10: The brick wall lattice with the assignment of Pauli matrices.

Derivation of an Euclidean action of a three-dimensional $(x,y,t)$, periodic
in all directions, system is very similar to the previous $(1+1)$-dimensional
example.

The two kinetic (i.e. $\sigma_{X}\sigma_{X}$ and $\sigma_{Y}\sigma_{Y}$) terms
in (43) are represented by the same six- or eight-spin couplings between the
adjacent time slices plus the appropriate phase, which naturally generalizes
the $(1+1)$-dimensional phase in the last term of (42) to three Euclidean
dimensions.

On the other hand diagonal, in the Hamiltonian form, potential terms are
represented by the standard Ising-like, ferromagnetic couplings along the
$y$-direction. They are located on the shorter edges of bricks at each time
slice. Hence, they are staggered in accord with the ($t$-independent) $x-y$
parity, $\zeta_{xy}=(-1)^{x+y}$, of a site originating given $Z$-edge in the
potential term. The final action reads

$\displaystyle S_{3D}$
$\displaystyle=\beta_{t}\sum_{x,y,t}O^{(6)}_{x,y,t}+\beta_{s}\sum_{\begin{subarray}{c}x,y,t,\\\
\zeta_{xy}=1\end{subarray}}O^{(2)}_{x,y,t}$ (44)
$\displaystyle+\frac{i\pi}{2}\sum_{\begin{subarray}{c}x,y,t,\\\
\zeta_{xy}=-1\end{subarray}}O^{(7)}_{x,y,t},$

with the phase operator $O^{(7)}_{x,y,t}$ being the direct generalization of
above $O^{(7)}_{x,t}$ to three dimensions and similarly for other couplings:

$\displaystyle
O^{(7)}_{x,y,t}=\frac{1}{2^{3}}(s_{x,y,t}+s_{x+1,y,t})(1+s_{x-1,y,t}s_{x-1,y,t+1})(1-s_{x,y,t}s_{x,y,t+1})(1-s_{x+1,y,t}s_{x+1,y,t+1}),$
(45a) $\displaystyle
O^{(6)}_{x,y,t}=\frac{1-2p}{8}(1+s_{x-1,y,t}s_{x-1,y,t+1})(1-s_{x,y,t}s_{x,y,t+1})(1-s_{x+1,y,t}s_{x+1,y,t+1})+\frac{p}{2}(1-s_{x,y,t}s_{x,y,t+1}),$
(45b) $\displaystyle O^{(2)}_{x,y,t}=-s_{x,y,t}s_{x,y+1,t}.$ (45c)

The action (44) describes then a three-dimensional Ising-like system. Together
with the corresponding constraints (still to be implemented) it would provide
an equivalent, Euclidean representation of a single, quantum Majorana spin on
a two-dimensional spatial lattice.

Even without the constraints the system is still interesting per se. Its
thermodynamics, the phase diagram, order parameters are unknown at the moment
and could be studied with standard methods of statistical physics. Such
studies would also provide, among other things, some information about the
constraints themselves.

The Boltzmann factor associated with (44) is not positive. However the origin
of its phases is now conceptually simple. Below we look how severe is the sign
problem in these unconstrained Euclidean models.

###  V.3 The sign problem

The standard (and practically only) method to deal with non-positive weights
is the reweighting [61, 62]. Instead of potentially negative Boltzmann factor
$\rho=\exp{(-S)}$, one uses as a Monte Carlo (MC) weight its absolute value
$\rho_{A}=|\rho|$, correcting at the same time all observables for this bias.

Whether such an approach is practical can be readily judged from the average
value of a sign of the exact Boltzmann factor

$\langle\mathrm{sign}\rangle\equiv\left\langle\frac{\rho}{\rho_{A}}\right\rangle_{A}=\frac{\mathcal{Z}}{\mathcal{Z}_{A}}$
(46)

averaged over the modulus $\rho_{A}$. If this average is close to $0$, the
method fails. If the contrary is true, say for some intermediate volumes, one
may expect to obtain meaningful estimates.

We have calculated analytically above average for both $(1+1)$\- and
$(2+1)$-dimensional models by employing the transfer matrix technique for a
range of small volumes. It is seen below that the sign problem is not very
severe in this case. Consequently, MC studies remain a viable approach to
explore these systems in detail.

Figure 11: Exact results for the average sign $\langle\mathrm{sign}\rangle$
for a range of two dimensional volumes $V$ and for the penalty parameter $p=2$
(left column) or $p=8$ (right column). Plots are presented for $\lambda$
values (from bottom to top) $0.1$, $0.25$, $0.5$, $0.75$, $1$, $1.25$, $1.5$,
$1.75$ and $2$.

####  V.3.1 $(1+1)$-dimensions

Partition functions $\mathcal{Z}$ and $\mathcal{Z}_{A}$ were calculated
exactly by summing Boltzmann factors $\exp{(-S_{2D})}$ and
$|\exp{(-S_{2D})}|$, as defined in Eqns (42). In Fig.11 the average sign is
shown for a range of two dimensional volumes and various couplings $\beta_{t}$
and $\beta_{s}$. The results are displayed as a function of a time step,
$\epsilon=\exp{(-\beta_{t})}$, and parameterized by different couplings
$\lambda=\frac{\beta_{s}}{\epsilon}$ in the Hamiltonian (39). Second column
displays analogous results for larger penalty parameter $p$.

The sign problem seems manageable for a sizeable part of the parameter space.
It vanishes entirely for $\epsilon\rightarrow 0$.

Increasing the penalty parameter $p$ also helps, since then some undesired
transitions vanish faster with $\epsilon$.

Both of these features show up also in our three-dimensional system. They can
be readily understood and used for our advantage, as discussed below.

(a)

(b)

Figure 12: Average sign as a function of the $\epsilon$ parameter in the
three-dimensional case with volume $V=4\times 4\times 3$ for the penalty
parameter (a) $p=2$ and (b) $p=8$. The range of $\lambda$ parameter is
explicitly given for both cases.

####  V.3.2 $(2+1)$-dimensions

For the three-dimensional Euclidean system (44) of volume $V=L_{x}L_{y}L_{t}$
a brute-force summation of all $2^{V}$ terms becomes already a challenge.
Still it was possible to obtain the value of $\langle\mathrm{sign}\rangle$ for
$V=4\times 4\times 3$, as shown in Fig. 12. It was done by constructing two
subsequent transfer matrices in the $y$ direction.

Again, as in the $(1+1)$-dimensions, the phase is harmless for small
$\epsilon$. This feature improves dramatically with increasing the penalty
parameter.

In addition, for $L_{t}=2$ no phase was observed in all cases. That is
$\langle\mathrm{sign}\rangle=1$ for all values of parameters and for all
studied dimensions.

####  V.3.3 The sign problem - summary

All the regularities observed above can be readily understood and generalized
for arbitrary sizes of lattices, providing at the same time some guidelines
for other, similar systems.

Consider first the case $L_{t}=2$. The partition function

$\mathcal{Z}^{(2)}=\sum_{s,s^{\prime}}e^{-\mathsf{L}(s,s^{\prime})}e^{-\mathsf{L}(s^{\prime},s)}$
(47)

is the sum over two-composite states of spins at the two time slices. The non-
zero phase can occur only if $s$ and $s^{\prime}$ differ by a double flip.
However in this case the phases of $e^{-\mathsf{L}(s,s^{\prime})}$ and
$e^{-\mathsf{L}(s^{\prime},s)}$ cancel and the result is positive for each
pair of configurations, as found above.

On the other hand, already for $L_{t}=3$ there are three states in the game

$\mathcal{Z}^{(3)}=\sum_{s,s^{\prime},s^{\prime\prime}}e^{-\mathsf{L}(s,s^{\prime\prime})-\mathsf{L}(s^{\prime\prime},s^{\prime})-\mathsf{L}(s^{\prime},s)}.$
(48)

Hence a single double-flip, e.g. in $s\rightarrow s^{\prime}$, can be balanced
by two subsequent single-flips in $s^{\prime}\rightarrow s^{\prime\prime}$ and
$s^{\prime\prime}\rightarrow s$ transitions. Since a phase may occur only in
the double flip transition $s\rightarrow s^{\prime}$, this particular
contribution may be negative and would give $\langle\mathrm{sign}\rangle<1$.

Consequently, the single flip transitions provide an undesired background
which indirectly causes negative signs of Boltzmann factors, hence the sign
problem.

However such transitions vanish for $\epsilon\rightarrow 0$ having a weight of
the higher order in $\epsilon$ by construction. This is clearly confirmed by
our calculations, cf. Figs 11 and 12, and explains why sign problem vanishes
at small $\epsilon$. These figures were obtained by using Mathematica [63].

Moreover, by increasing the penalty parameter $p$ we can force the “bad
transitions” to vanish faster. Indeed this is also confirmed by our results
for $p=8$ in both dimensions. This suggests that the sign problem could be
significantly reduced by setting $p=\infty$, which amounts to introducing a
constraint in the Euclidean system 222Not to be confused with plaquette
constraints required for bosonization..

We remark that even under this constraint there exist “Euclidean histories”
with negative sign. As they involve a number of spin flips growing with the
system size, one may hope that they do not lead to significant difficulties.

Obviously all these scenarios should be further studied quantitatively.

##  VI Conclusions and outlook

We have presented a bosonization method generalizing the idea from [4], valid
for lattices of arbitrary coordination number and with arbitrary number of
Majorana modes per lattice site. In the previous works only systems with even
coordination numbers and one pair of fermionic creation/annihilation operators
per lattice site were considered. The new approach extends the construction in
several ways. First, for lattices with vertices of even degree we may include
multiple fermionic states per site. We illustrate this by bosonizing the
Hubbard model. Second, we allow for lattices with odd coordination numbers.
Then there is an odd number of Majorana fermions per site. In particular
systems with one Majorana per site may be bosonized.

It is possible to treat also systems for which the coordination number is not
congruent modulo two to the number of Majoranas per site. Strictly speaking in
this case we do not bosonize the original fermionic system but rather one
augumented by some spurious fermionic degrees of freedom. Nevertheless,
operators corresponding on the bosonic side to these modes may be clearly
identified and decoupled. Even in the case of very regular lattices such trick
is needed in presence of a boundary. We emphasize that this is a feature of
our bosonization method, not of fermionic systems per se.

If these two numbers are not congruent modulo two, it involves spurious
fermionic states, which nevertheless can be identified and eliminated. Another
potential source of interest in this construction is that it provides new
analytically tractable examples of spin systems featuring edge modes.

One question which remains unanswered is whether our construction may be
dualized to some higher gauge theory. For systems with one fermion (and hence
two Majoranas) per lattice site such picture of bosonization has been obtained
in [2, 3].

Concerning the Euclidean formulation, our main conclusion is that in spite of
somewhat unusual time evolution, generated by simultaneous double-flips, a
local Euclidean action for an unconstrained system was derived. It contains at
least $6$-spin interactions and is highly asymmetric between space and time,
in contrast to the standard Ising model. To our knowledge, this system has not
been studied. Now, it can be readily explored with standard statistical
methods.

Our generic action is complex. It was found that the resulting sign problem is
manageable on small lattices where our fully analytical approach is available.

The next logical step now is to study the problem for larger, although
intermediate, sizes and see whether the popular reweighting methods allow
meaningful measurements of observables, extrapolation to larger volumes and
extraction of scaling limits.

Moreover, it is conceivable that introduced here methods could be extended to
implement the spin constraints avoiding the standard non-positive Legendre
transformation.

We intend to further study some of these questions with the aid of
quantitative Monte Carlo approach.

###### Acknowledgements.

The authors acknowledge W. Brzezicki for suggesting the decagonal lattice as a
three-dimensional example with odd degree vertices. AB acknowledges G. Ortiz
and Z. Nussinov for pointed out the references [50, 51, 52] and for discussion
about them. This work is supported in part by the NCN grant:
UMO-2016/21/B/ST2/01492. BR was also supported by the MNS donation for PhD
students and young scientists N17/MNS/000040.

## References

  * Jordan and Wigner [1928] P. Jordan and E. Wigner, Über das paulische äquivalenzverbot, Z. Phys. 47, 631 (1928).
  * Chen  _et al._ [2018a] Y.-A. Chen, A. Kapustin, and D. Radičević, Exact bosonization in two spatial dimensions and a new class of lattice gauge theories, Annals of Physics 393, 234 (2018a).
  * Chen and Kapustin [2019] Y.-A. Chen and A. Kapustin, Bosonization in three spatial dimensions and a 2-form gauge theory, Phys. Rev. B 100, 245127 (2019).
  * Wosiek [1982] J. Wosiek, A local representation for fermions on a lattice, Acta Phys. Polon. B 13, 543 (1982).
  * Bochniak and Ruba [2020] A. Bochniak and B. Ruba, Bosonization based on Clifford algebras and its gauge theoretic interpretation, JHEP 12, 118, arXiv:2003.06905 [math-ph] .
  * Burgess  _et al._ [1994] C. Burgess, C. Lütken, and F. Quevedo, Bosonization in higher dimensions, Physics Letters B 336, 18 (1994).
  * Bravyi and Kitaev [2002] S. B. Bravyi and A. Y. Kitaev, Fermionic quantum computation, Annals of Physics 298, 210 (2002).
  * Ball [2005] R. C. Ball, Fermions without fermion fields, Phys. Rev. Lett. 95, 176407 (2005).
  * Verstraete and Cirac [2005] F. Verstraete and J. I. Cirac, Mapping local hamiltonians of fermions to local hamiltonians of spins, Journal of Statistical Mechanics: Theory and Experiment 2005, P09012 (2005).
  * Fradkin [1989] E. Fradkin, Jordan-wigner transformation for quantum-spin systems in two dimensions and fractional statistics, Phys. Rev. Lett. 63, 322 (1989).
  * Karch and Tong [2016] A. Karch and D. Tong, Particle-vortex duality from 3d bosonization, Phys. Rev. X 6, 031043 (2016).
  * Zohar and Cirac [2018] E. Zohar and J. I. Cirac, Eliminating fermionic matter fields in lattice gauge theories, Phys. Rev. B 98, 075119 (2018).
  * Chen [2020] Y.-A. Chen, Exact bosonization in arbitrary dimensions, Phys. Rev. Research 2, 033527 (2020).
  * Son [2015] D. T. Son, Is the composite fermion a dirac particle?, Phys. Rev. X 5, 031027 (2015).
  * Chen  _et al._ [2018b] J.-Y. Chen, J. H. Son, C. Wang, and S. Raghu, Exact boson-fermion duality on a 3d euclidean lattice, Phys. Rev. Lett. 120, 016602 (2018b).
  * Mattis and Lieb [1965] D. C. Mattis and E. H. Lieb, Exact solution of a many‐fermion system and its associated boson field, Journal of Mathematical Physics 6, 304 (1965), https://doi.org/10.1063/1.1704281 .
  * Witten [1984] E. Witten, Non-abelian bosonization in two dimensions, Commun. Math. Phys. 92, 455 (1984).
  * Cobanera  _et al._ [2011] E. Cobanera, G. Ortiz, and Z. Nussinov, The bond-algebraic approach to dualities, Advances in Physics 60, 679 (2011), https://doi.org/10.1080/00018732.2011.619814 .
  * Kramers and Wannier [1941] H. A. Kramers and G. H. Wannier, Statistics of the two-dimensional ferromagnet. part i, Phys. Rev. 60, 252 (1941).
  * Franz [1971] W. Franz, Duality in generalized ising models and phase transitions without local order parameter, Journal of Mathematical Physics 12 (1971).
  * Karch  _et al._ [2019] A. Karch, D. Tong, and C. Turner, A Web of 2d Dualities: ${\bf Z}_{2}$ Gauge Fields and Arf Invariants, SciPost Phys. 7, 7 (2019).
  * Seiberg  _et al._ [2016] N. Seiberg, T. Senthil, C. Wang, and E. Witten, A duality web in 2+1 dimensions and condensed matter physics, Annals of Physics 374, 395 (2016).
  * Senthil  _et al._ [2019] T. Senthil, D. T. Son, C. Wang, and C. Xu, Duality between $(2+1)d$ quantum critical points, Physics Reports 827, 1 (2019), duality between $(2+1)d$ quantum critical points.
  * Schultz  _et al._ [1964] T. D. Schultz, D. C. Mattis, and E. H. Lieb, Two-dimensional ising model as a soluble problem of many fermions, Rev. Mod. Phys. 36, 856 (1964).
  * Mandal and Surendran [2009] S. Mandal and N. Surendran, Exactly solvable kitaev model in three dimensions, Phys. Rev. B 79, 024426 (2009).
  * Condella and DeTar [2000] J. Condella and C. DeTar, Potts flux tube model at nonzero chemical potential, Phys. Rev. D 61, 074023 (2000).
  * Delgado  _et al._ [2014] Y. Delgado, C. Gattringer, and A. Schmidt, Solving the sign problem of two flavor scalar electrodynamics at finite chemical potential, PoS LATTICE 2013, 147 (2014).
  * Gattringer  _et al._ [2015] C. Gattringer, T. Kloiber, and V. Sazonov, Solving the sign problems of the massless lattice schwinger model with a dual formulation, Nuclear Physics B 897, 732 (2015).
  * Kitaev [2003] A. Kitaev, Fault-tolerant quantum computation by anyons, Annals of Physics 303, 2 (2003).
  * Kitaev and Laumann [2009] A. Kitaev and C. Laumann, Topological phases and quantum computation (2009), arXiv:0904.2771 [cond-mat.mes-hall] .
  * Gaiotto and Kapustin [2016] D. Gaiotto and A. Kapustin, Spin TQFTs and fermionic phases of matter, Int. J. Mod. Phys. A 31, 1645044 (2016), arXiv:1505.05856 [cond-mat.str-el] .
  * Mazaheri  _et al._ [2015] T. Mazaheri, G. Ortiz, Z. Nussinov, and A. Seidel, Zero modes, bosonization, and topological quantum order: The laughlin state in second quantization, Phys. Rev. B 91, 085115 (2015).
  * Kapustin and Thorngren [2017] A. Kapustin and R. Thorngren, Fermionic SPT phases in higher dimensions and bosonization, JHEP 10, 080, arXiv:1701.08264 [cond-mat.str-el] .
  * Cappelli  _et al._ [2017] A. Cappelli, E. Randellini, and J. Sisti, Three-dimensional topological insulators and bosonization, J. High Energ. Phys. 135.
  * Pozo  _et al._ [2021] O. Pozo, P. Rao, C. Chen, and I. Sodemann, Anatomy of $\mathbb{Z}_{2}$ fluxes in anyon fermi liquids and bose condensates, Phys. Rev. B 103, 035145 (2021).
  * Szczerba [1985] A. M. Szczerba, Spins and fermions on arbitrary lattices, Commun. Math. Phys. 98, 513 (1985).
  * Bochniak  _et al._ [2020] A. Bochniak, B. Ruba, J. Wosiek, and A. Wyrzykowski, Constraints of kinematic bosonization in two and higher dimensions, Phys. Rev. D 102, 114502 (2020), arXiv:2004.00988 [hep-lat] .
  * Majorana [1937] E. Majorana, Teoria simmetrica dell’elettrone e del positrone, Nuovo Cim 14, 171 (1937).
  * Note [1] More precisely, complex structure compatible with the bilinear form determining the canonical anticommutation relations.
  * Schrieffer [1964] J. R. Schrieffer,  _Theory of Superconductivity_ (W. A. Benjamin, 1964).
  * Bardeen  _et al._ [1957] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Theory of superconductivity, Phys. Rev. 108, 1175 (1957).
  * Abrikosov [1957] A. A. Abrikosov, On the Magnetic properties of superconductors of the second group, Sov. Phys. JETP 5, 1174 (1957).
  * Wilczek [2009] F. Wilczek, Majorana returns, Nature Phys , 614–618 (2009).
  * Read and Green [2000] N. Read and D. Green, Paired states of fermions in two dimensions with breaking of parity and time-reversal symmetries and the fractional quantum hall effect, Phys. Rev. B 61, 10267 (2000).
  * Kopnin and Salomaa [1991] N. B. Kopnin and M. M. Salomaa, Mutual friction in superfluid ${}^{3}\mathrm{He}$: Effects of bound states in the vortex core, Phys. Rev. B 44, 9667 (1991).
  * Moore and Read [1991] G. Moore and N. Read, Nonabelions in the fractional quantum hall effect, Nuclear Physics B 360, 362 (1991).
  * Kitaev [2001a] A. Y. Kitaev, Unpaired majorana fermions in quantum wires, Physics-Uspekhi 44, 131 (2001a).
  * Sarma  _et al._ [2015] S. Sarma, M. Freedman, and C. Nayak, Majorana zero modes and topological quantum computation, npj Quantum Inf 1, 15001 (2015).
  * Huang  _et al._ [2021] H.-L. Huang, M. Narożniak, F. Liang, Y. Zhao, A. D. Castellano, M. Gong, Y. Wu, S. Wang, J. Lin, Y. Xu, H. Deng, H. Rong, J. P. Dowling, C.-Z. Peng, T. Byrnes, X. Zhu, and J.-W. Pan, Emulating quantum teleportation of a majorana zero mode qubit, Phys. Rev. Lett. 126, 090502 (2021).
  * Nussinov and Ortiz [2009] Z. Nussinov and G. Ortiz, Bond algebras and exact solvability of hamiltonians: Spin $s=\frac{1}{2}$ multilayer systems, Phys. Rev. B 79, 214440 (2009).
  * Yao  _et al._ [2009] H. Yao, S.-C. Zhang, and S. A. Kivelson, Algebraic spin liquid in an exactly solvable spin model, Phys. Rev. Lett. 102, 217202 (2009).
  * Wu  _et al._ [2009] C. Wu, D. Arovas, and H.-H. Hung, $\Gamma$-matrix generalization of the kitaev model, Phys. Rev. B 79, 134427 (2009).
  * Wilczek [1982] F. Wilczek, Magnetic flux, angular momentum, and statistics, Phys. Rev. Lett. 48, 1144 (1982).
  * Kitaev [2006] A. Kitaev, Anyons in an exactly solved model and beyond, Annals Phys. 321, 2 (2006), arXiv:cond-mat/0506438 .
  * Hubbard [1963] J. Hubbard, Electron correlations in narrow energy bands, Proc. R. Soc. Lond. A 276, 238–257 (1963).
  * Nielsen and Ninomiya [1981] H. Nielsen and M. Ninomiya, Absence of neutrinos on a lattice: (i). proof by homotopy theory, Nuclear Physics B 185, 20 (1981).
  * Friedan [1982] D. Friedan, A proof of the nielsen-ninomiya theorem, Communications in Mathematical Physics 85, 481 (1982).
  * Kitaev [2001b] A. Y. Kitaev, Unpaired majorana fermions in quantum wires, Physics-Uspekhi 44, 131 (2001b).
  * Fradkin and Susskind [1978] E. Fradkin and L. Susskind, Order and disorder in gauge systems and magnets, Phys. Rev. D 17, 2637 (1978).
  * Kogut [1979] J. B. Kogut, An introduction to lattice gauge theory and spin systems, Rev. Mod. Phys. 51, 659 (1979).
  * Landau and Binder [2014] D. P. Landau and K. Binder,  _A Guide to Monte Carlo Simulations in Statistical Physics_ (Cambridge University Press, 2014).
  * Nakamura  _et al._ [1992] T. Nakamura, N. Hatano, and H. Nishimori, Reweighting method for quantum monte carlo simulations with the negative-sign problem, Journal of the Physical Society of Japan 61, 3494 (1992).
  * [63] W. R. Inc., System Modeler, Version 12.1, champaign, IL, 2020\.
  * Note [2] Not to be confused with plaquette constraints required for bosonization.