# XXZ spin chain at half-filling subject to local dissipation

Oscar Bouverot-Dupuis  Université Paris Saclay, CNRS, LPTMS, 91405, Orsay,
France  IPhT, CNRS, CEA, Université Paris Saclay, 91191 Gif-sur-Yvette, France
Saptarshi Majumdar  Université Paris Saclay, CNRS, LPTMS, 91405, Orsay, France
IPhT, CNRS, CEA, Université Paris Saclay, 91191 Gif-sur-Yvette, France
Alberto Rosso  Université Paris Saclay, CNRS, LPTMS, 91405, Orsay, France
Laura Foini  IPhT, CNRS, CEA, Université Paris Saclay, 91191 Gif-sur-Yvette,
France

###### Abstract

We study an XXZ spin chain at zero magnetization coupled to a collection of
local harmonic baths at zero temperature. We map this system on a (1+1)D
effective field theory using bosonization, where the effect of the bath is
taken care of in an exact manner. We provide analytical and numerical evidence
of the existence of two phases at zero temperature: a Luttinger liquid (LL)
and an antiferromagnetic phase (AFM), separated by a
Berezinsky–Kosterlitz–Thouless (BKT) phase transition. While the bath is
responsible for the LL-AFM transition for subohmic baths, the LL-AFM
transition for superohmic baths is due to the interactions within the spin
chain.

##  I Introduction

The problem of a quantum-mechanical system interacting with an environment is
ubiquitous in physics. These interactions can originate from a surrounding
bath, an external driving force, an optical or solid lattice, or else.
Regardless, the resulting behavior of the system coupled to its environment
can drastically differ from that of the isolated system. At one end of the
spectrum, if the environment evolves on time scales much larger than that of
the system, it is said to be frozen. In the case where the environment settles
in a random frozen state, it can be accurately modeled by quenched disorder.
This is the case of the celebrated Anderson model [1] where a particle gets
localized due to the presence of a random potential, leading to vanishing
transport properties. The other end of the spectrum corresponds to
environments with very rapid dynamics that can safely be considered Markovian
and for which the Lindblad master equation has proven very useful [2]. In
between these two time scales are slowly varying baths for which a successful
quantum formulation was first established by Caldeira, Leggett et al. [3, 4].
They studied a single degree of freedom (spin or particle) coupled to a bath,
such as the spin-boson model [5] or a quantum Brownian particle in a periodic
potential [6, 7]. In these models, the system is coupled to a bath of phonons
that is simple enough to be traced over and recover effective dynamics for the
system of interest.

More recently, interest has started shifting towards many-body versions of
these models. Going from one particle to many opens the door to a new zoology
of phenomena, but some single-particle effects can also persist. For instance,
adding interactions to the Anderson model is conjectured to preserve
localization at finite temperature through the so-called Many-Body
Localization (MBL) [8, 9, 10], i.e. the fact that some isolated many-body
systems can fail to reach equilibrium and retain a memory of their initial
conditions.

In this work, we study the well-known XXZ spin chain coupled to local baths of
harmonic oscillators (see fig. 1). To a certain extent, this model can be
viewed as a many-body extension in (1+1)D of the spin-boson model in (0+1)D
[5, 11], or an extension of the most studied many-body localization setup, an
XXZ spin chain in a random magnetic field [8], where the fields are now
replaced by local baths. The effect of the bath on its corresponding spin can
be captured by an exponent $s$ (see the next section for a more formal
definition) which allows to model a large variety of baths. For a spin chain
at finite magnetization and any $s$, it has been shown that the bath induces
fractional excitations that create a dissipative phase presenting signatures
of localization [12, 13, 14, 15]. At zero magnetization, the only known
results are for $s=1$, i.e. an ohmic bath, for which a superfluid to insulator
transition has been proposed in [16, 15]. In this article, we broaden the
scope of these studies by focusing on the zero magnetization case and allowing
for a generic exponent $s$. Our results show the existence of a
Berezinsky–Kosterlitz–Thouless (BKT) phase transition between a Luttinger
liquid (LL) and an antiferromagnetic phase (AFM). The exact location of the
transition point depends on the bath exponent $s$ and is shifted compared to
the non-dissipative spin chain for $s<1$.

The manuscript is organized as follows: section II introduces the model and
maps it to a bosonic effective field theory. Section III then presents a clear
overview of the main results in terms of this effective field-theoretic
description. These results were derived using thorough numerical and
analytical approaches. On the analytical side, a perturbative renormalization
group study is presented in section IV and is complemented by a variational
approach in section V. These analytical predictions are then tested for
several observables against an exact numerical simulation of the field theory
in section VI. Finally, a brief discussion of the results and concluding
remarks are made in section VII. Unless specified otherwise (as in section V),
the zero temperature ($\beta\to\infty$) and thermodynamic ($L\to\infty$)
limits will always be understood.

##  II Model

Figure 1: Schematic representation of the model. An XXZ spin chain with
parameters $J_{\rm z},J_{\rm xy}$ has its spins coupled to independent and
identical collections of harmonic oscillators. The spin-bath interaction is
assumed to be fully described by the parameters $\alpha$ and $s$.

The XXZ spin chain is a 1D periodic chain of N spins, total length $L$, and
lattice spacing $a=\frac{L}{N}$. It evolves according to the Hamiltonian

$H_{\rm S}=\sum_{j=1}^{N}J_{\rm z}S_{j}^{z}S_{j+1}^{z}-J_{\rm
xy}\left(S_{j}^{x}S_{j+1}^{x}+S_{j}^{y}S_{j+1}^{y}\right),$ (1)

where $S^{\mu}_{i}$ are the spin-1/2 operators. Each spin is then coupled to a
set of quantum harmonic oscillators by

$H_{\rm SB}=\sum_{j=1}^{N}S_{j}^{z}\sum_{k}\lambda_{k}X_{k,j},$ (2)

and the oscillators’ Hamiltonian is

$H_{\rm
B}=\sum_{j=1}^{N}\sum_{k}\frac{P_{k,j}^{2}}{2m_{k}}+\frac{1}{2}m_{k}\Omega_{k}^{2}X_{k,j}^{2}.$
(3)

Hence the total Hamiltonian of the dissipative system $H=H_{\rm S}+H_{\rm
SB}+H_{\rm B}$. Note that the baths are local, i.e. each bath has its
independent spin that it acts upon and the baths are not spatially correlated.
Thus, its effect can be thought of as an instance of annealed disorder
$h_{j}(t)=\sum_{k}\lambda_{k}X_{k,j}$, acting locally through $H_{\rm
SB}=\sum_{j=1}^{N}h_{j}(t)S_{j}^{z}$. The quenched limit of this problem
$h_{j}(t)=h_{j}$ has already been studied in [17, 18] where a zero temperature
delocalization-localization transition was found.

A common tool to study such many-body systems is the Lindblad master equation
[2] which makes use of a Markovian approximation of the bath’s dynamics.
Nevertheless, in this work, we will not use it and resort to an approach
similar to that of Caldeira and Leggett [5], which takes the bath into account
in an exact manner. Following the idea of [5], the bath characteristics are
completely encoded in the spectral function

$J(\Omega)=\frac{\pi}{2}\sum_{k}\frac{\lambda_{k}^{2}}{\Omega_{k}m_{k}}\delta(\Omega-\Omega_{k}),$
(4)

which we assume to have the low-energy behavior

$J(\Omega)=\alpha\frac{\pi\Omega_{D}^{1-s}}{\Gamma(1+s)}\Omega^{s}\text{ for
}\Omega\in[0,\Omega_{D}],$ (5)

where $s$ is the bath exponent and $\alpha$ measures the strength of the
coupling to the bath. Following the classification from Caldeira and Leggett,
we call $s=1$ an ohmic bath, while $0<s<1$ is a subohmic bath and $s>1$ is a
superohmic bath. The rest of this section will be dedicated to deriving an
effective action for our system using the bosonization technique.

###  II.1 Bosonization

To implement the bosonization procedure described in [19, 20, 12], one must
first map the spin chain to a one-dimensional spinless fermionic system using
the Jordan–Wigner transformation. With the usual notations
$c_{j},c_{j}^{\dagger}$ for the ladder operators and
$n_{j}=c_{j}^{\dagger}c_{j}$, the XXZ spin chain Hamiltonian becomes $H_{\rm
S}=H_{\rm xy}+H_{\rm z}$, where

$\displaystyle H_{\rm xy}$ $\displaystyle=-\frac{J_{\rm
xy}}{2}\sum_{j=1}^{N}\left(c^{\dagger}_{j+1}c_{j}+c^{\dagger}_{j}c_{j+1}\right),$
(6) $\displaystyle H_{\rm z}$ $\displaystyle=\sum_{j=1}^{N}J_{\rm
z}\left(n_{j}-\frac{1}{2}\right)\left(n_{j+1}-\frac{1}{2}\right).$ (7)

The hopping Hamiltonian $H_{\rm xy}$ can be diagonalized in Fourier space as
$H_{\rm xy}=\sum_{k}\varepsilon_{k}c^{\dagger}_{k}c_{k}$ with
$\varepsilon_{k}=-J_{\rm xy}\cos(ka)$. In the rest of the work, we will
consider the spin chain to be at zero magnetization, which corresponds to a
half-filled fermionic chain and a Fermi momentum $k_{F}=\frac{\pi}{2a}$ that
is commensurate with the lattice spacing $a$. The application of the
bosonization technique then requires linearizing the spectrum
$\varepsilon_{k}$ around $k_{F}$ which leads to a field-theoretic description
of the Hamiltonian $H_{\rm S}$ with the action $S_{\rm S}=S_{\rm LL}+S_{g}$,
given by

$\displaystyle S_{\rm LL}$ $\displaystyle=$
$\displaystyle\int\frac{dxd\tau}{2\pi
K}\left[u\left(\partial_{x}\phi(x,\tau)\right)^{2}+\frac{1}{u}\left(\partial_{\tau}\phi(x,\tau)\right)^{2}\right]\;\;\;$
(8) $\displaystyle S_{\rm g}$ $\displaystyle=$
$\displaystyle-\frac{gu}{2\pi^{2}a^{2}}\int dxd\tau\cos[4\phi(x,\tau)],$ (9)

where $x\in[0,L]$ denotes the spatial coordinate, $\tau\in[0,\beta]$ is the
imaginary time coordinate with $\beta$ the inverse temperature of the system,
and $S_{\rm LL}$ describes a Luttinger liquid (LL). The parameters $u,K$ are
the so-called Luttinger parameters and, with $g$, are related to the
microscopic parameters. The Bethe ansatz gives their exact expressions (e.g.
$K^{-1}_{\rm Bethe}=\frac{2}{\pi}\arccos(-J_{\rm z}/J_{\rm xy})$) while the
bosonization yields approximate expressions valid in the $J_{\rm z}\ll J_{\rm
xy}$ limit (e.g. $K^{-1}_{\rm bosonization}=\sqrt{1+\frac{4J_{\rm z}}{\pi
J_{\rm xy}}}$). On top of this description of the isolated spin chain, one
needs to add the effect of the bath captured by $H_{\rm SB}+H_{\rm B}$. This
requires the bosonized expression of $S_{j}^{z}$

$S^{z}_{j}=-\frac{a}{\pi}\nabla\phi(x_{j})+\frac{(-1)^{j}}{\pi}\cos[2\phi(x_{j})],\>x_{j}=ja.$
(10)

Using this expression to write a path integral representation of $H_{\rm
B}+H_{\rm SB}$, one realises that the bath degrees of freedom are quadratic
and can therefore be integrated out. This leads to the following additional
term in the action

$\displaystyle S_{\alpha}=$ $\displaystyle-\frac{a}{2\pi^{2}}\int dxd\tau
d\tau^{\prime}$
$\displaystyle\times\left[\partial_{x}\phi(x,\tau)-\frac{(-1)^{x/a}}{a}\cos[2\phi(x,\tau)]\right]\mathcal{D}(\tau,\tau^{\prime})$
$\displaystyle\times\left[\partial_{x}\phi(x,\tau^{\prime})-\frac{(-1)^{x/a}}{a}\cos[2\phi(x,\tau^{\prime})]\right],$
(11)

where one can show using the spectral function $J(\Omega)$ that the
dissipative kernel is
$\mathcal{D}(\tau,\tau^{\prime})\sim\alpha\tau_{c}^{s-1}|\tau-\tau^{\prime}|^{-1-s}$
for $\tau>\tau_{c}=1/\Omega_{D}$ (see appendix A for the detailed
computation). For the sake of simplicity, we relate the imaginary-time and
space cutoffs as $a=u\tau_{c}$ which doesn’t change the underlying physics.
Eq. II.1 can be further simplified by dropping the rapidly fluctuating terms
$(-1)^{x/a}\partial_{x}\phi(x,\tau)\cos[2\phi(x,\tau^{\prime})]$, and the
cross gradient term
$\partial_{x}\phi(x,\tau)\partial_{x}\phi(x,\tau^{\prime})$ which is
irrelevant by power counting. The total bosonized action $S=S_{\rm
LL}+S_{g}+S_{\alpha}$ is thus

$\displaystyle S[\phi]=$ $\displaystyle\int dxd\tau\frac{1}{2\pi
K}\left[u\left(\partial_{x}\phi(x,\tau)\right)^{2}+\frac{1}{u}\left(\partial_{\tau}\phi(x,\tau)\right)^{2}\right]-\frac{g}{2\pi^{2}a\tau_{c}}\int
dxd\tau\cos[4\phi(x,\tau)]$ (12)
$\displaystyle-\frac{\alpha}{2\pi^{2}a\tau_{c}^{1-s}}\underset{|\tau-\tau^{\prime}|>\tau_{c}}{\int
dxd\tau
d\tau^{\prime}}\frac{\cos[2\phi(x,\tau)]\cos[2\phi(x,\tau^{\prime})]}{|\tau-\tau^{\prime}|^{1+s}}.$

This effective action is that of a Luttinger liquid (LL) with two types of
interactions: a local sine-Gordon interaction controlled by $g$, and a long-
range interaction controlled by $\alpha$ which is commonly encountered in
generalized XY models [21, 22] or dissipative quantum systems [23, 24]. Note
that the long-range interaction involving
$\cos[2\phi(x,\tau)]\cos[2\phi(x,\tau^{\prime})]$ can be rewritten as the sum
of $\cos[2(\phi(x,\tau)+\phi(x,\tau^{\prime}))]$ and
$\cos[2(\phi(x,\tau)-\phi(x,\tau^{\prime}))]$. In the case of an
incommensurate XXZ spin chain [12, 13, 14], only the latter term remains and
the sine-Gordon term vanishes. It has been shown that this term gives rise to
fractional excitations $|\omega_{n}|^{s}$ in the spectrum of the
incommensurate spin chain. However, from our variational analysis in section
V, we will demonstrate that this fractional term is absent from our
commensurate model and is replaced by a gap.

##  III Main results

From the bosonized action in Eq. 12, the analytical and numerical tools used
in the following sections infer the zero temperature ($\beta\to\infty$) and
thermodynamic ($L\to\infty$) phase diagram depicted in fig. 2. A Luttinger
liquid (LL) and an antiferromagnet (AFM) are separated by a transition that
depends on the bath exponent $s$.

  * •

For a superohmic bath ($s>1$), a BKT transition 111Here and in the following,
we call BKT transition a transition whose RG equations are of the BKT type.
occurs at $K_{c}=1/2$ (in the limit of infinitesimal $\alpha,g$) and is driven
by the coupling $g$ which comes from the internal interactions of the spin
chain. The spin chain is an LL for $K\geq K_{c}$ while it becomes an AFM at
$K<K_{c}$.

  * •

For a subohmic bath ($s<1$), the system-bath interaction, denoted by the
coupling $\alpha$, shifts the BKT transition to $K_{c}=1-s/2$ (in the limit of
infinitesimal $\alpha,g$). This increased value of $K_{c}$ extends the area of
the AFM.

In both cases, the critical point $K_{c}$ increases with $\alpha$ and $g$. The
two phases are characterized as follows. On the one hand, the LL is a critical
phase exhibiting quasi-long range order as seen by the power-law decay of
spin-spin correlation functions [19]. On the other hand, the AFM is a gapped
(or massive) phase exhibiting long-range order. Although suggested by [16], we
do not find any fractional excitations in the AFM phase, neither from the
analytical study, nor from the exact numerical simulations. This AFM is of the
same nature as the AFM in the sine-Gordon model. The transition admits the
order parameter $\langle\cos(2\phi)\rangle$ which is related to the
antiferromagnetic spin density wave as $\langle
S_{j}^{z}\rangle=\frac{(-1)^{j}}{\pi}\langle\cos(2\phi)\rangle$. We show that
$\langle\cos(2\phi)\rangle$ vanishes with the gap of the AFM phase, thus
making it an order parameter of the AFM-LL transition.

Figure 2: Phase diagram in the $(K,s)$ plane as described by the RG analysis
in sec. IV and the variational method in sec. V for infinitesimal value of
$\alpha$ and $g$. The golden-colored phase is a gapless LL, and the blue-
colored phase is a gapped AFM phase. Along the dashed line, the phase
transition is controlled by the coupling $\alpha$, while along the dotted
line, it is governed by $g$.

##  IV Perturbative renormalization group approach

Looking at the bosonized action in Eq. (12), we expect that, at large
distances and times, the model presents at least two phases: a Luttinger
liquid (LL) phase where both couplings $\alpha$ and $g$ are irrelevant, and a
strongly interacting one where the couplings are relevant. To capture the
precise location of the departure from the LL, we implement a perturbative
renormalization group (RG) analysis. It turns out it is enough to compute the
RG equations up to $\mathcal{O}(\alpha,g^{2},\alpha g)$ to grasp the
interesting physics at play. The perturbative RG analysis was done using the
operator product expansion formalism [16, 26] to respect the real-space sharp
cutoff $a$ appearing in the action. The detailed computation can be found in
appendix B and leads to the following RG equations where the dependence on the
renormalization time $l$ has been made explicit,

$\displaystyle\frac{d}{dl}$
$\displaystyle\frac{u(l)}{K(l)}=\frac{g(l)^{2}u(l)}{\pi^{2}},$ (13)
$\displaystyle\frac{d}{dl}$
$\displaystyle\frac{1}{u(l)K(l)}=\frac{2\alpha(l)}{\pi
u(l)}+\frac{g(l)^{2}}{\pi^{2}u(l)},$ (14) $\displaystyle\frac{d}{dl}$
$\displaystyle g(l)=(2-4K(l))g(l)+\alpha(l),$ (15) $\displaystyle\frac{d}{dl}$
$\displaystyle\alpha(l)=(2-s-2K(l))\alpha(l)+\frac{g(l)\alpha(l)}{\pi},$ (16)

which agree with [16, 12]. Notice how the presence of $\alpha$ in the RG
equation for $g$ implies that the bath will generate the coupling $g$ even if
its microscopic (bare) value is 0. The symmetric process, i.e. $g$ generates
$\alpha$, is not possible as $g$ can only enhance a non zero $\alpha$ as seen
in Eq. 16. When considering the effect of the bath (i.e. $\alpha\neq 0$), this
means that there are not two distinct phase transitions at $K_{c}=1/2$ and
$K_{c}=1-s/2$ as seen from the scaling dimensions of $\alpha$ and $g$, but
rather a unique transition at $K_{c}=\max(1-s/2,1/2)$. The corresponding phase
diagram is depicted in fig. 2 and shows the existence of a LL for $K>K_{c}$,
with $K_{c}=1-s/2$ for a subohmic bath ($s<1$) and $K_{c}=1/2$ for a
superohmic bath ($s>1$). Note that the long-range cosine potential is
irrelevant by power counting for $s>2$. In this work, we are interested in
understanding the effect of the bath on the system, thus we constrain the
value of $s\in[0,2]$ in fig. 2.

###  IV.1 Superohmic bath ($s>1$)

For a superohmic bath ($s>1$), near the transition point $K_{c}=1/2$, the
coupling $\alpha$ has a scaling dimension $2-s-2K_{c}=1-s<0$ which is strongly
irrelevant. It can therefore be safely ignored and the RG equations become

$\displaystyle\frac{d}{dl}$ $\displaystyle u=0,$ (17)
$\displaystyle\frac{d}{dl}$ $\displaystyle\frac{1}{K}=\frac{g^{2}}{\pi^{2}},$
(18) $\displaystyle\frac{d}{dl}$ $\displaystyle g=(2-4K)g,$ (19)

which are the one-loop RG equations of the sine-Gordon model [27, 19]. The
associated RG flow is shown in fig. 3 and is known to belong to the BKT
universality class.

Figure 3: (left) Superohmic transition: RG flow of the couplings $g$ and $K$
according to Eqs. 18,19. We distinguish two phases: a Luttinger liquid (LL)
depicted in gold, and an antiferromagnetic phase (AFM) in blue. The critical
point is at $K_{c}=1/2$. (right) Subohmic transition (here $s=0.5$): RG flow
of the couplings $\alpha$ and $K$ according to Eqs. 23,24. We distinguish two
phases: a Luttinger liquid (LL) in gold, and an antiferromagnetic phase (AFM)
in blue. The critical point is at $K_{c}=1-s/2=0.75$.

###  IV.2 Ohmic bath ($s=1$)

The ohmic bath ($s=1$) is the most studied type of bath in the literature [12,
16, 15]. In this setting, the scaling dimensions of $\alpha$ and $g$ both
vanish at $K_{c}=1/2$, which means the transition is driven by both couplings
simultaneously. The RG equations Eqs. 13,14,15,16 cannot be simplified any
further (as is the case for superohmic or subohmic baths). This leads to an RG
flow similar to that of the BKT transition, the main difference being that it
is governed by two couplings ($\alpha$ and $g$) instead of one.

###  IV.3 Subohmic bath ($s<1$)

For a subohmic bath ($s<1$), the coupling $g$ is strongly irrelevant near the
transition point $K_{c}=1-s/2$ since its scaling dimension is
$2-4K_{c}=2s-2<0$. The RG equations can therefore be simplified by discarding
$g$ to give

$\displaystyle\frac{d}{dl}$ $\displaystyle\frac{u}{K}=0,$ (20)
$\displaystyle\frac{d}{dl}$ $\displaystyle\frac{1}{uK}=\frac{2\alpha}{\pi u},$
(21) $\displaystyle\frac{d}{dl}$ $\displaystyle\alpha=(2-s-2K)\alpha.$ (22)

The LL action appears to be renormalized along the imaginary-time direction
($\frac{d}{dl}\frac{1}{uK}\neq 0$) but not along the space direction
($\frac{d}{dl}\frac{u}{K}=0$)222Pushing the RG equations to higher orders
shows that $\frac{d}{dl}\frac{u}{K}=\mathcal{O}(\alpha^{2})$. This is to be
contrasted with the incommensurate spin chain [12, 13] for which
$\frac{d}{dl}\frac{u}{K}=0$ is non-perturbative. This is due to the fact that
the system-bath interaction in the incommensurate spin chain action is
invariant under the statistical tilt symmetry $\phi(x,\tau)\to\phi(x,\tau)+\mu
x$, where $\mu$ is a constant.. This anisotropy is due to the long-range
interaction in Eq. 12 which only spans the imaginary-time direction. Since the
ratio $\frac{u}{K}$ is constant, Eqs. 21,22 can be written in terms of $K$ and
$\sqrt{\alpha}$ only

$\displaystyle\frac{d}{dl}$
$\displaystyle\frac{1}{K}=\frac{\sqrt{\alpha}^{2}}{\pi},$ (23)
$\displaystyle\frac{d}{dl}$
$\displaystyle\sqrt{\alpha}=(1-s/2-K)\sqrt{\alpha},$ (24)

which are exactly the same type of equations as Eqs. 18,19 describing the
sine-Gordon model. This shows that this phase transition is also of the BKT
type and doesn’t belong to a new universality class as suggested by [29]. The
RG flow associated to these equations is depicted in fig. 3.

##  V Variational method

While the perturbative RG is tailored to capture the end of the Luttinger
liquid (LL), it fails at describing correctly the physics of the
antiferromagnetic phase (AFM). This is why we resort to a variational method à
la Feynman [30, 19] to capture the essential properties of the bulk of the
phases. This technique aims at finding an approximate action $S_{\rm var}$
that is very similar to $S$ while being simple enough to allow for analytical
computations. We consider a quadratic ansatz

$S_{\rm
var}=\frac{1}{2}\int\frac{dq}{2\pi}\frac{d\omega_{n}}{2\pi}\phi^{\star}(q,\omega_{n})G^{-1}_{\rm
var}(q,\omega_{n})\phi(q,\omega_{n}),$ (25)

where $G^{-1}_{\rm var}(q,\omega_{n})=1/G_{\rm var}(q,\omega_{n})$ is to be
determined and $\phi(q,\omega_{n})$ is the Fourier transform of $\phi(x,\tau)$
given by
$\phi(x,\tau)=\int\frac{dq}{2\pi}\frac{d\omega_{n}}{2\pi}\phi(q,\omega_{n})e^{i\left(qx-\omega_{n}\tau\right)}$.
The distance from the original action to the variational action is then
defined through the variational free energy $F_{\rm var}=-T\log Z_{\rm
var}+T\langle S-S_{\rm var}\rangle_{\rm var}$ ($\langle\ldots\rangle_{\rm
var}$ stands for the average with respect to $S_{\rm var}$). Indeed, the true
free energy of the original field theory $F=-T\log Z=T\log Z_{\rm
var}-T\log\langle e^{-(S-S_{var})}\rangle_{\rm var}$, is always upper bounded
by $F_{\rm var}$ due to the convexity of the exponential. This defines the
best variational action as that which minimizes $F_{\rm var}$. The goal is
thus to solve the vanishing gradient condition $\frac{\delta F_{\rm
var}}{\delta G^{-1}_{\rm var}(q,\omega_{n})}=0$. It is possible to give an
explicit expression for $F_{\rm var}$ which in turn yields the vanishing
gradient equation

$\displaystyle G^{-1}_{\rm var}(q,\omega_{n})=$ $\displaystyle\frac{1}{\pi
K}\left[uq^{2}+\frac{\omega_{n}^{2}}{u}\right]+\frac{8g}{\pi^{2}a\tau_{c}}\exp\left[-\frac{2}{\pi^{2}}\int
dq^{\prime}d\omega_{n}^{\prime}G_{\rm
var}(q^{\prime},\omega_{n}^{\prime})\right]$ (26)
$\displaystyle+\frac{2\alpha}{\pi^{2}a\tau_{c}^{1-s}}\int_{\tau_{c}}^{\infty}\frac{d\tau}{\tau^{1+s}}\sum_{\varepsilon=\pm
1}(1+\varepsilon\cos\omega_{n}\tau)\exp\left[-\frac{1}{\pi^{2}}\int
dq^{\prime}d\omega_{n}^{\prime}(1+\varepsilon\cos\omega_{n}^{\prime}\tau)G_{\rm
var}(q^{\prime},\omega_{n}^{\prime})\right].$

###  V.1 Phase diagram at infinitesimal coupling

In the absence of interactions, i.e. $\alpha=g=0$, Eq. 26 reduces to
$G^{-1}_{\rm LL}(q,\omega_{n})=\frac{1}{\pi
K}\left[\frac{\omega_{n}^{2}}{u}+uq^{2}\right]$ which is just the propagator
of the Luttinger liquid. When reintroducing the interaction terms, we expect a
gap to appear, like in the sine-Gordon model [19], and maybe some fractional
excitations, as in [16], so we postulate the following ansatz for the
propagator in the AFM

$\displaystyle G^{-1}_{\rm AFM}(q,\omega_{n})=$ $\displaystyle\frac{1}{\pi
K}\left[uq^{2}+\frac{\omega_{n}^{2}}{u}+\nu|\omega_{n}|^{s}+\frac{\Delta^{2}}{u}\right],$
(27)

where the parameters $\Delta$ and $\nu$ are to be determined. This ansatz is
expected to be valid in the $\alpha,g\to 0$ limit where one can neglect the
renormalization of the Luttinger parameters $u,K$ into $u_{r},K_{r}$ (more on
this in the next subsection). The first observation we make is that $\nu=0$.
This can be understood from the fact that the $|\omega_{n}|^{s}$ excitations
can only come from the bath-dependent $\alpha$ term in Eq. 27, and the
$\varepsilon=1$ contribution exactly cancels out that of $\varepsilon=-1$. It
is worth noticing that unlike [16], where a similar action with an ohmic bath
was studied, this means that we do not find any fractional excitation
$|\omega_{n}|^{s}$ in our system for any $s$. We now move on to the
determination of the gap $\Delta$. Setting $q=0$, $\omega_{n}=0$ in Eq. 26
yields an equation for $\Delta$ that can be written as (see appendix C for the
detailed computation)

$\displaystyle\frac{(\Delta\tau_{c})^{2-4K}}{K}=$ $\displaystyle
g\frac{8}{\pi}+\alpha\frac{4e^{2K\gamma_{E}}}{\pi}\frac{(\Delta\tau_{c})^{s-2K}-1}{2K-s},$
(28)

where $\tau_{c}$ is the imaginary-time UV cutoff, and $\gamma_{E}$ is Euler’s
gamma constant. This equation is valid in the limit of small $\Delta\tau_{c}$
which we expect to happen when $K\sim K_{c}$. To arrive at the phase diagram
shown in fig. 2, we write Eq. 28 as
$\Delta^{2-4K}=a_{1}+a_{2}\frac{\Delta^{s-2K}}{2K-s}$ with $a_{1},a_{2}$ two
constants independent of $\Delta$. As we approach the transition $K\to
K_{c}^{-}$, the gap $\Delta$ is expected to vanish, so we retain only the
leading terms in the previous gap equation. This suggests to distinguish three
cases depending on the sign of $2K_{c}-s$.

Figure 4: Numerical evaluation of observables for the commensurate dissipative
spin chain. The values of the parameters are taken to be $s=0.5$, $K=1$,
$g=0$, $u=1$, and $\beta$ is scaled as $L$. For the left and middle column,
red, green, and blue colors denote $L=128$, $L=256$ and $L=320$ respectively.
All the quantities have been averaged over 10000-20000 configurations and the
value of time-step $dt$ is $0.5$. (top) and (bottom) rows correspond to
$\alpha=1$ (LL) and $\alpha=5$ (AFM), respectively. (left) Susceptibility
$\chi$ as a function $q$. For $\alpha=1$, it stays finite and constant whereas
for $\alpha=5$, it vanishes at small $q$. (middle) Behavior of $C(\omega_{n})$
as a function of $\omega_{n}$. For $\alpha=1$, $\omega_{n}C(\omega_{n})$
saturates to a constant $K_{r}/2=0.477$; for $\alpha=5$, $C(\omega_{n})$
becomes a constant $K_{r}/(2\Delta_{r})=0.155$. (right) Behavior of the order
parameter $\langle\cos\left[2\left(\phi-\phi_{\rm CoM}\right)\right]\rangle$
as a function of system size $L$. For $\alpha=1$, the order parameter decays
as $L^{-K_{r}}$ with $K_{r}=0.933$, whereas for $\alpha=5$, it increases and
saturates to a constant as $c_{1}-c_{2}/L$ with $c_{1}=0.679$ and
$c_{2}=0.894$.

  * •

$s-2K_{c}>0$: The gap equation reduces to $\Delta=a_{1}^{\frac{1}{2-4K}}$
which tells us that the critical point is $K_{c}=1/2$ with $K<K_{c}$
corresponding to the gapped phase, while $K>K_{c}$ is gapless. We also
understand that this gapped phase solution is valid only if $s>2K_{c}=1$, i.e.
the bath is superohmic ($s>1$).

  * •

$s-2K_{c}<0$: Retaining only the leading terms yields $\Delta\sim
a_{2}^{\frac{1}{2-s-2K}}$, signifying that $K_{c}=1-\frac{s}{2}$ and the
gapped phase is for $K<K_{c}$. In this regime $s<2K_{c}=2-s$, which means the
bath is subohmic ($s<1$).

  * •

$s-2K_{c}=0$: For $K$ close to $K_{c}$, one expands $a^{2K-s}=1+(2K-s)\log a$
since $2K-s\ll 1$, which shows that to leading order
$\Delta^{2-4K}\propto-\log(\Delta\tau_{c})$. This yields $K_{c}=1/2$ and thus
$s=1$ which is an ohmic bath. This matches with the two previous cases but the
gap is now discontinuous at the transition. This is a known artifact of the
variational method [16, 19].

In the spin chain picture, this gapped phase turns out to be an AFM. Indeed,
by averaging Eq. 10 over all field configurations, one arrives at

$\langle
S^{z}_{j}\rangle=\frac{(-1)^{\frac{x_{j}}{a}}}{\pi}\langle\cos[2\phi(x_{j})]\rangle,$
(29)

which denotes the existence of an AFM phase with amplitude
$\frac{1}{\pi}\langle\cos(2\phi)\rangle$. Using the variational propagator
$G^{-1}_{\rm AFM}$ for the gapped phase, one then proves that
$\langle\cos(2\phi)\rangle=\left(\Delta\tau_{c}\right)^{K}$, which can be
obtained from the $\Delta\tau_{c}\ll 1$ limit of Eq. 74. This implies that
$\langle\cos(2\phi)\rangle$ can be used as an order parameter for this
transition as this quantity vanishes in the LL and is finite in the AFM.
Moreover, looking back at the gap expressions derived near the transition, it
appears that all derivatives of the order parameter vanish at the transition.
This corroborates the scenario of a BKT phase transition which is known to be
of infinite order. Physically, this order parameter is associated with the
spontaneous symmetry breaking of the discrete shift symmetry
$\phi(x,\tau)\to\phi(x,\tau)+\frac{n\pi}{2}$ of the bosonized field. Indeed,
we expect $\langle\cos(2\phi)\rangle=\langle\cos(2\phi+\pi)\rangle=0$ when the
symmetry holds, which shows that the AFM is the symmetry breaking phase.

The fact this order parameter is associated to an AFM phase is not surprising
as in the case of a dissipative incommensurate spin chain, the ordered phase
is a spin density wave of wavelength $\pi/k_{F}$ [13]. Putting
$k_{F}=\frac{\pi}{2a}$, one recovers an AFM spin density wave with a $2a$
wavelength. All these results confirm the phase diagram obtained from the RG
(see fig. 2). There are two phases separated by $K_{c}=\max(1-s/2,1/2)$; for
$K>K_{c}$, the system remains an LL, while for $K<K_{c}$ it becomes an AFM.

###  V.2 Phase diagram at finite coupling

In the previous subsection, we showed how the simple ansatz $G^{-1}_{\rm
AFM}(q,\omega_{n})=\frac{1}{\pi
K}\left[uq^{2}+\frac{\omega_{n}^{2}}{u}+\frac{\Delta^{2}}{u}\right]$ was
enough to capture the correct location of the phase transition when
$\alpha,g\to 0$. However, from the RG results we expect that, starting at any
value of $K>K_{c}$ and with $g$ fixed, the system remains in the LL phase for
$0<\alpha<\alpha_{c}$ and switches to an AFM for $\alpha>\alpha_{c}$. At
$\alpha=\alpha_{c}$, the system should be a LL with renormalized Luttinger
parameter $K_{r}=K_{c}$. The same analysis is of course valid if tuning $g$
while keeping $\alpha$ constant. One could therefore wonder if a more general
ansatz $G^{-1}_{\rm var}(q,\omega_{n})=\frac{1}{\pi
K_{r}}\left[u_{r}q^{2}+\frac{\omega_{n}^{2}}{u_{r}}+\frac{\Delta_{r}^{2}}{u_{r}}\right]$,
where $u_{r}$, $K_{r}$, $\Delta_{r}$ are fitting parameters to be determined,
could capture the renormalization of the Luttinger parameters as predicted by
the RG. It appears that this ansatz works well deep in the phase for finite
$\alpha$. However, it fails to recover the lowest order RG equations close to
the transition (see appendix D). In the next section, we will nonetheless use
the respective renormalized quantities, denoted with a subscript $r$ (e.g. the
renormalized value of $K$ is $K_{r}$ and so on), as fitting parameters for our
numerical analysis done at finite $g,\alpha$.

##  VI Numerical Results

In this section, we calculate different observables, both analytically and
numerically, to observe the signature of the phase transition and physically
characterize the ordered phase (AFM). For our numerical analysis, we simulate
the Langevin dynamical equation for the field $\phi$ associated to the
equilibrium probability distribution $P_{\rm
eq.}[\phi]=e^{-S\left[\phi\right]}$. The Langevin equation is
$\frac{d\phi_{ij}(t)}{dt}=-\frac{\delta
S\left[\phi_{ij}(t)\right]}{\delta\phi_{ij}(t)}+\eta_{ij}(t)$, where $i,j$ are
the discretized indices for imaginary-time $\tau$ and space $x$ respectively,
$t$ denotes the Langevin time (alternatively, the simulation time) and
$\eta(t)$ is Gaussian white noise with $\langle\eta_{ij}(t)\rangle=0$ and
$\langle\eta_{ij}(t)\eta_{i^{\prime}j^{\prime}}(t^{\prime})\rangle=2\delta_{i,i^{\prime}}\delta_{j,j^{\prime}}\delta(t-t^{\prime})$.
Note that the noise $\eta$ is used only for thermalizing the Langevin equation
and is not to be confused with the noise coming from the dissipative bath,
which has already been taken into account in the action $S[\phi]$. Using this
numerical technique, we simulate configurations for fixed values of $s,K$ and
for different values of $\alpha$. For each set of parameters, we extract field
configurations of different sizes $L\times\beta$. We scale $L$ as $\beta$,
keeping the BKT dynamic scaling $z=1$ in mind; and on these configurations, we
calculate a number of observables to characterize the two phases. In the
following, we show the results for $s=0.5$ (subohmic bath), $K=1$, $g=0$, and
$u=1$. Results for $s=1$ (ohmic bath) and $s=1.5$ (superohmic bath) can be
found in appendix G.

###  VI.1 Relevant observables

The first quantity that we calculate is the static susceptibility of the spin
chain $\chi=\lim_{q\to 0}(q/\pi)^{2}G\left(q,\omega_{n}=0\right)$, also known
as the compressibility of the field. Using the LL propagator $G_{\rm LL}=\pi
K_{r}\left[u_{r}q^{2}+\frac{\omega_{n}^{2}}{u_{r}}\right]^{-1}$ derived in
section V, it appears that the static susceptibility is finite in the LL phase
and is given by $\chi_{\rm LL}=K_{r}/(u_{r}\pi)$. On the other hand, using
$G_{\rm AFM}=\pi
K_{r}\left[u_{r}q^{2}+\frac{\omega_{n}^{2}}{u_{r}}+\frac{\Delta_{r}^{2}}{u_{r}}\right]^{-1}$
shows that $\chi_{\rm AFM}$ vanishes as
$\frac{K_{r}u_{r}}{\Delta_{r}^{2}}q^{2}$ in the gapped phase and for $q\to 0$.
These analytical predictions are then compared to the data from the Langevin
simulation. Numerically, the Green’s function $G\left(q,\omega_{n}\right)$ is
obtained by computing the correlation function
$\langle\phi(q,\omega_{n})\phi(-q,-\omega_{n})\rangle=\langle|\phi(q,\omega_{n})|^{2}\rangle$,
as the fields are real and bosonic. The associated numerical results are
depicted in fig. 4, left and support our analytical predictions. In the LL
($\alpha=1$), the quantity $\lim_{q\to
0}(q^{2}/\pi)G\left(q,\omega_{n}=0\right)=\pi\chi$ is equal to $K_{r}/u_{r}=1$
(top row). On the other hand, in the AFM ($\alpha=5$), the $q\to 0$ limit
vanishes (bottom row).

Figure 5: Numerical phase diagram obtained for $K=1,u=1,g=0,s=0.5$ by varying
$\alpha$. (left) $K_{r}$ in the LL phase as a function of $\alpha$, extracted
from the order parameter $\langle\cos\left[2\left(\phi-\phi_{\rm
CoM}\right)\right]\rangle$ (purple circles) and from $C(\omega_{n})$ (blue
triangles). $K_{r}$ decreases to $K_{c}=0.75$ as $\alpha$ is increased,
signalling a BKT transition. (right) As $\alpha$ is increased, the gap (red
dots) first becomes finite for $\alpha=3$ and then increases.

The next quantity that we calculate is given by $C(\omega_{n})=\frac{1}{\pi
L}\sum_{q=-\infty}^{\infty}G(q,\omega_{n})$. In the thermodynamic limit, the
sum over $q$ can be replaced with the integral
$\frac{1}{2\pi}\int_{-\infty}^{\infty}dq$. In this limit, using our
variational action from section V we expect that for a small-$\omega_{n}$
limit,

$C(\omega_{n}\to 0)=\begin{cases}\frac{K_{r}}{2\omega_{n}}\ &\mbox{LL,}\\\
\frac{K_{r}}{2\Delta_{r}}\ &\mbox{AFM.}\end{cases}$ (30)

Comparing this to the numerical results in fig. 4, middle, we see that indeed
for small $\alpha$ (LL phase), $\omega_{n}C(\omega_{n})$ goes to a constant
for $\omega_{n}\to 0$; whereas $C(\omega_{n})$ itself saturates to a constant
for a higher value of $\alpha$ (AFM phase). This confirms the existence of a
gap in the low-energy spectrum of the dissipative phase.

Finally, we compute the order parameter
$\langle\cos\left[2\left(\phi-\phi_{\rm CoM}\right)\right]\rangle$, where
$\phi_{\rm CoM}=\frac{1}{\beta L}\sum_{x,\tau}\phi(x,\tau)$ is the center of
mass (CoM) of the configuration. The field has been offset by its CoM as to
suppress the contribution $G(0,0)$ to the order parameter which would
otherwise diverge in the gapless LL phase. In appendix F, we analytically show
that in the LL phase (with $\beta=L$), the order parameter vanishes as
$L^{-K_{r}}$, whereas it increases and saturates to a constant in the AFM as
$c_{1}-c_{2}/L$, where $c_{1}$ and $c_{2}$ are both positive constants.
Numerical results in fig. 4, right, corroborate this scenario. This confirms
our prediction that the dissipative bath can induce a quantum phase transition
on a 1D XXZ spin chain by spontaneously breaking the discrete shift symmetry
$\phi\to\phi+\frac{n\pi}{2},n\in\mathcal{Z}$, and the ordered phase is a
gapped AFM phase.

###  VI.2 Microscopic parameters

One can numerically extract relevant microscopic parameters in both phases
from fitting the previous observables. In the LL phase, $K_{r}$ can be
extracted from the fittings of the order parameter and $C(\omega_{n})$. In the
gapped phase, the extraction of the parameters is slightly more tricky as, at
lowest order, the quantitative behavior of the phase is regulated by the gap
term. By computing $\lim_{q\to 0}\chi(q)/q^{2}$, one can nonetheless extract
the value of the inverse of the gap $u_{r}K_{r}/\Delta_{r}^{2}$. The behaviors
of these parameters are given in fig. 5. From the plots, we see that $K_{r}$
starts decreasing as $\alpha$ is increased and approaches $K_{c}=0.75$ at the
transition, after which $\Delta_{r}$ becomes finite and increases as $\alpha$
is increased. This tells us that for $K=1$, $u=1$, $g=0$ and $s=0.5$, the
transition happens around $\alpha_{c}\in(2,3)$.

##  VII Discussions and conclusions

In this work, we investigated the zero temperature properties of a XXZ spin
chain at zero magnetisation coupled to local baths of phonons _à la_ Caldeira
and Leggett. Using the bosonization procedure, this system was mapped on to an
effective field theory (with bosonic field $\phi$) which could thus be tackled
using powerful analytical and numerical field-theoretic tools. We confirmed
the existence of a BKT phase transition between a Luttinger liquid (LL)
exhibiting quasi-long-range order and an antiferromagnetic (AFM) phase with
long-range order. The variational method predicts that the AFM does not
present any fractional excitations as suggested by [16]. The exact location of
the phase transition depends on the type of bath studied through the bath
exponent $s$. For superohmic baths ($s>1$), dissipation does not affect the
transition which remains that of the isolated spin chain to leading order.
However, for subohmic baths ($s<1$), the transition is shifted and the AFM
eats into the LL. It appears that this transition is associated with the
spontaneous symmetry breaking of the order parameter
$\langle\cos(2\phi)\rangle$. This quantity identifies with the amplitude of
the antiferromagnetic spin density wave
$S^{z}_{j}\sim(-1)^{j}\langle\cos(2\phi)\rangle$. As with a regular BKT
transition, all derivatives of this order parameter vanish at the transition,
signaling an infinite order phase transition.

While the variational method used in section V successfully captured the main
properties of both phases, it makes use of a rather crude approximation
amounting to replacing the system’s highly non-linear action by a quadratic
one. It is a well-established fact that such an approximate action fails at
capturing topological excitations of the system such as solitons or instantons
[19, 31]. For a generic AFM with a vanishing linear conductivity, such
excitations might lead to the restoring of non-linear terms [32, 33]. For our
model there are no instantons with finite action in the zero temperature limit
(see chapter 3.2 of [31] for a proof of this statement), and the system is
truly locked in one of the minima of the potential. Nonetheless, one might
worry about the effect of solitons, for instance on the transport properties,
as for the pure sine-Gordon model. We leave for further investigation the
understanding of this point. Moreover, the study of the same model at finite
temperature, where instantons connecting degenerate minima play an important
role, seems an interesting path to explore.

Another direction we wish to explore is that of the competition between bath-
induced localization, as studied in this work and [12, 13], and disorder-
induced Anderson localization. One such example can be found in [34] where a
quantum phase transition was found between an Anderson localized phase and a
Zeno localized phase in a one-dimensional non-interacting system. Finally,
valuable insights could probably be gained by mapping this system to a
Coulomb-like gas using the general idea of the sine-Gordon to the 2D Coulomb
gas mapping [35].

###### Acknowledgements.

Acknowledgments: This work was granted access to the HPC resources of IDRIS
under the allocation 2022-(AD011013581R1) made by GENCI. We thank Thierry
Giamarchi, Thibaud Maimbourg, and Manon Michel for useful discussions. We are
also grateful to Edmond Orignac for pointing out the role of topological
solutions in our system.

##  Appendix A Dissipative kernel

When integrating out the bath degrees of freedom, one generates the
dissipative kernel
$\mathcal{D}(\tau,\tau^{\prime})=\sum\limits_{k}\mathcal{D}_{k}(\tau,\tau^{\prime})$,
where
$\mathcal{D}_{k}^{-1}(\tau,\tau^{\prime})=\frac{m_{k}}{\lambda_{k}^{2}}\delta(\tau-\tau^{\prime})(\Omega_{k}^{2}-\partial_{\tau^{\prime}}^{2})$.
To recover $\mathcal{D}(\tau,\tau^{\prime})$, one must therefore invert
$\mathcal{D}^{-1}_{k}(\tau,\tau^{\prime})$ which amounts to finding
$\mathcal{D}_{k}(\tau,\tau^{\prime})$ such that

$\int
d\tau^{\prime}\mathcal{D}_{k}(\tau,\tau^{\prime})\mathcal{D}^{-1}_{k}(\tau,\tau^{\prime})=\delta(\tau-\tau^{\prime\prime}).$
(31)

Using the following Fourier and inverse transform conventions

$\displaystyle\mathcal{D}_{k}(\omega,\omega^{\prime})$ $\displaystyle=\int
d\tau
d\tau^{\prime}e^{i\omega\tau}\mathcal{D}_{k}(\tau,\tau^{\prime})e^{-i\omega^{\prime}\tau^{\prime}},$
(32) $\displaystyle\mathcal{D}_{k}(\tau,\tau^{\prime})$
$\displaystyle=\int\frac{d\omega}{2\pi}\frac{d\omega^{\prime}}{2\pi}e^{-i\omega\tau}\mathcal{D}_{k}(\omega,\omega^{\prime})e^{i\omega^{\prime}\tau^{\prime}},$
(33)

Eq. 31 can be written in Fourier space as

$\int
d\omega^{\prime}\mathcal{D}_{k}(\omega,\omega^{\prime})\mathcal{D}_{k}^{-1}(\omega^{\prime},\omega^{\prime\prime})=4\pi^{2}\delta(\omega-\omega^{\prime\prime}).$
(34)

Let’s apply this to the kernel $\mathcal{D}_{k}^{-1}(\tau,\tau^{\prime})$.
According to Eq. 32, its Fourier transform is
$\mathcal{D}^{-1}_{k}(\omega,\omega^{\prime})=\frac{m_{k}}{\lambda_{k}^{2}}2\pi\delta(\omega-\omega^{\prime})(\Omega_{k}^{2}+{\omega^{\prime}}^{2})$.
From this expression and Eq. 34, it is clear that
$\mathcal{D}_{k}(\omega,\omega^{\prime})=\frac{\lambda_{k}^{2}}{m_{k}}\frac{2\pi\delta(\omega-\omega^{\prime})}{\Omega_{k}^{2}+{\omega^{\prime}}^{2}}$.
This implies, along with Eq. 33, that

$\displaystyle\mathcal{D}_{k}(\tau,\tau^{\prime})$
$\displaystyle=\frac{\lambda_{k}^{2}}{m_{k}}\int\frac{d\omega}{2\pi}\frac{e^{-i\omega(\tau-\tau^{\prime})}}{\Omega_{k}^{2}+\omega^{2}}=\frac{\lambda_{k}^{2}}{2m_{k}\Omega_{k}}e^{-\Omega_{k}|\tau-\tau^{\prime}|}.$
(35)

Performing the sum over $k$ then yields

$\displaystyle\mathcal{D}(\tau,\tau^{\prime})$ $\displaystyle=\int
d\Omega\,e^{-\Omega|\tau-\tau^{\prime}|}\sum_{k}\frac{\lambda_{k}^{2}}{2m_{k}\Omega_{k}}\delta(\Omega-\Omega_{k}),$
(36)

where we recognize the spectral function from Eq. 4 as
$J(\Omega)=\frac{\pi}{2}\sum_{k}\frac{\lambda_{k}^{2}}{m_{k}\Omega_{k}}\delta(\Omega-\Omega_{k})=\alpha\frac{\pi\Omega_{D}^{1-s}}{\Gamma(1+s)}\Omega^{s}\text{
for }\Omega\in[0,\Omega_{D}]$. This leads to

$\displaystyle\mathcal{D}(\tau,\tau^{\prime})$
$\displaystyle=\alpha\frac{\Omega_{D}^{1-s}}{\Gamma(1+s)}\int_{0}^{\Omega_{D}}d\Omega\,\Omega^{s}e^{-\Omega|\tau-\tau^{\prime}|}$
$\displaystyle=\frac{\alpha\Omega_{D}^{1-s}}{|\tau-\tau^{\prime}|^{1+s}}\frac{\int_{0}^{\Omega_{D}|\tau-\tau^{\prime}|}dx\,x^{s}e^{-x}}{\Gamma(1+s)}.$
(37)

For $|\tau-\tau^{\prime}|\gg\tau_{c}=1/\Omega_{D}$, one can approximate the
integral as $\int_{0}^{\infty}dxx^{s}e^{-x}=\Gamma(1+s)$. Thus

$\mathcal{D}(\tau,\tau^{\prime})=\frac{\alpha\tau_{c}^{s-1}}{|\tau-\tau^{\prime}|^{1+s}},$
(38)

which is the expression given in the main text.

##  Appendix B Perturbative RG analysis

The goal of this appendix is to compute the RG equations Eqs. 13,14,15,16 by
means of the operator product expansion (OPE).

###  B.1 Derivation of the useful OPEs

The operator product expansion (OPE) is a series expansion of a product of two
nearby fields in terms of local fields. This is done on normal ordered
operators to avoid any divergence not coming from the two operators being
pushed together. For our problem, we will need the following OPEs

$\displaystyle\mathrel{\mathop{\mathchar 58\relax}}e^{2ip\phi(r)}$
$\displaystyle\mathrel{\mathop{\mathchar 58\relax}}\mathrel{\mathop{\mathchar
58\relax}}e^{2ip\phi(r^{\prime})}\mathrel{\mathop{\mathchar
58\relax}}=a^{2p^{2}K}\mathrel{\mathop{\mathchar
58\relax}}e^{4ip\phi(r)}\mathrel{\mathop{\mathchar 58\relax}},$ (39)
$\displaystyle\mathrel{\mathop{\mathchar 58\relax}}e^{2ip\phi(r)}$
$\displaystyle\mathrel{\mathop{\mathchar 58\relax}}\mathrel{\mathop{\mathchar
58\relax}}e^{-2ip\phi(r^{\prime})}\mathrel{\mathop{\mathchar
58\relax}}=\frac{\mathrel{\mathop{\mathchar
58\relax}}e^{2ip(\phi(r)-\phi(r^{\prime}))}\mathrel{\mathop{\mathchar
58\relax}}}{|\delta r|^{2p^{2}K}}$
$\displaystyle=\frac{\mathrel{\mathop{\mathchar 58\relax}}1+2ip\delta
r\cdot\nabla\phi(r)-2p^{2}(\delta
r\cdot\nabla\phi(r))^{2}\mathrel{\mathop{\mathchar 58\relax}}}{|\delta
r|^{2p^{2}K}},$ (40)

which imply

$\displaystyle\mathrel{\mathop{\mathchar
58\relax}}\cos[2p\phi(r)]\mathrel{\mathop{\mathchar 58\relax}}$
$\displaystyle\mathrel{\mathop{\mathchar
58\relax}}\cos[2p\phi(r^{\prime})]\mathrel{\mathop{\mathchar
58\relax}}=\frac{a^{2p^{2}K}}{2}\mathrel{\mathop{\mathchar
58\relax}}\cos[4p\phi(r)]\mathrel{\mathop{\mathchar 58\relax}}$
$\displaystyle-p^{2}\frac{\mathrel{\mathop{\mathchar 58\relax}}(\delta
r\cdot\nabla\phi(r))^{2})\mathrel{\mathop{\mathchar 58\relax}}}{|\delta
r|^{2p^{2}K}}+...,$ (41) $\displaystyle\mathrel{\mathop{\mathchar
58\relax}}\cos[4\phi(r)]\mathrel{\mathop{\mathchar 58\relax}}$
$\displaystyle\mathrel{\mathop{\mathchar
58\relax}}\cos[2\phi(r^{\prime})]\mathrel{\mathop{\mathchar
58\relax}}=\frac{\mathrel{\mathop{\mathchar
58\relax}}\cos[2\phi(r)]\mathrel{\mathop{\mathchar 58\relax}}}{2a^{4K}}+...$
(42)

where $r=(x,u\tau)$, $r^{\prime}=(x^{\prime},u\tau^{\prime})$, $\delta
r=r^{\prime}-r$. These relations are easily proven using the identity
$\mathrel{\mathop{\mathchar 58\relax}}e^{i2p\phi(r)}\mathrel{\mathop{\mathchar
58\relax}}=\frac{e^{i2p\phi(r)}}{\langle
e^{i2p\phi(r)}\rangle}=\frac{e^{i2p\phi(r)}}{a^{p^{2}K}}$, where
$\langle\cdot\rangle$ is the average with respect to the Gaussian action
$S_{\rm LL}$.

###  B.2 Renormalization group using the OPE

The interacting part of the action in Eq. 12 can be rewritten in terms of
normal ordered fields as

$\displaystyle S_{\rm int}=$ $\displaystyle S_{g}+S_{\alpha}$ $\displaystyle=$
$\displaystyle-\frac{gu}{2\pi^{2}a^{2-4K}}\int
dxd\tau\mathrel{\mathop{\mathchar
58\relax}}\cos\left[4\phi(x,\tau)\right]\mathrel{\mathop{\mathchar 58\relax}}$
$\displaystyle-\frac{\alpha
u^{1-s}}{2\pi^{2}a^{2-s-2K}}\underset{|\tau-\tau^{\prime}|>\tau_{c}}{\int
dxd\tau d\tau^{\prime}}$ $\displaystyle\times\frac{\mathrel{\mathop{\mathchar
58\relax}}\cos\left[2\phi(x,\tau)\right]\mathrel{\mathop{\mathchar
58\relax}}\mathrel{\mathop{\mathchar
58\relax}}\cos\left[2\phi(x,\tau^{\prime})\right]\mathrel{\mathop{\mathchar
58\relax}}}{|\tau-\tau^{\prime}|^{1+s}}.$ (43)

One then expands the partition function up to order
$\mathcal{O}(\alpha,g^{2},\alpha g)$ such that

$\displaystyle Z$ $\displaystyle=\int\mathcal{D}\phi e^{-S_{\rm
LL}[\phi]-S_{g}[\phi]-S_{\alpha}[\phi]}$ $\displaystyle=Z_{\rm
LL}\Big{[}1-\langle S_{g}\rangle(a)-\langle S_{\alpha}\rangle(a)$
$\displaystyle\hskip 28.45274pt+\frac{\langle S_{g}^{2}\rangle(a)}{2}+\langle
S_{\alpha}S_{g}\rangle(a)\Big{]},$ (44)

where the cutoff dependence has been made explicit. We now perform a rescaling
of the lattice spacing $a\to a^{\prime}=a(1+dl)$ (so
$\tau_{c}\to\tau_{c}^{\prime}=\tau_{c}(1+dl)$) and ask how the couplings
should vary to preserve the partition function $Z$. $Z_{\rm LL}$ being an RG
fixed-point, we only need to consider the variations of the averaged
interacting terms. Keeping only the first order in $dl$, $\langle
S_{g}\rangle(a^{\prime})$ becomes

$\displaystyle\langle S_{g}\rangle(a^{\prime})$ $\displaystyle=\langle
S_{g}\rangle(a)\left[1+\frac{dg}{g}+(4K-2)dl\right].$ (45)

The next term $\langle S_{\alpha}\rangle(a^{\prime})$ requires a bit more
work. We start by splitting the imaginary-time integral as
$\underset{|\tau-\tau^{\prime}|>\tau_{c}}{\int[...]}-\underset{\tau_{c}^{\prime}>|\tau-\tau^{\prime}|>\tau_{c}}{\int[...]}$.
The first term is simply given by

$\displaystyle\langle
S_{\alpha}\rangle(a)\left[1+\frac{d\alpha}{\alpha}+(2K+s-2)dl\right],$ (46)

while the second term, which involves only fields at very close positions, can
be evaluated using the OPE B.1 and gives

$\displaystyle dl$
$\displaystyle\Big{[}\frac{\alpha(a)u}{2\pi^{2}a^{2-4K}}\int
dxd\tau\langle\mathrel{\mathop{\mathchar
58\relax}}\cos\left[4\phi(x,\tau)\right]\mathrel{\mathop{\mathchar
58\relax}}\rangle$ $\displaystyle-\frac{\alpha(a)}{\pi^{2}u}\int
dxd\tau\langle\mathrel{\mathop{\mathchar
58\relax}}(\partial_{\tau}\phi(x,\tau))^{2}\mathrel{\mathop{\mathchar
58\relax}}\rangle\Big{]}.$ (47)

The third term, namely $\frac{1}{2}\langle S_{g}^{2}\rangle(a^{\prime})$,
reads

$\displaystyle\frac{\langle
S_{g}^{2}\rangle(a^{\prime})}{2}=\frac{g^{2}a^{\prime 8K-4}}{32\pi^{4}}$
$\displaystyle\underset{|r-r^{\prime}|>a^{\prime}}{\int drdr^{\prime}}$ (48)
$\displaystyle\times\langle\mathrel{\mathop{\mathchar 58\relax}}\cos
4\phi(r)\mathrel{\mathop{\mathchar 58\relax}}\mathrel{\mathop{\mathchar
58\relax}}\cos 4\phi(r^{\prime})\mathrel{\mathop{\mathchar 58\relax}}\rangle.$

For this term, one again splits the limits of the integral into
$\int\underset{|r-r^{\prime}|>a}{drdr^{\prime}}[...]-\int\underset{a^{\prime}>|r-r^{\prime}|>a}{drdr^{\prime}}[...].$
The first part is easily tractable while the second requires the use of the
OPE B.1. This yields

$\displaystyle\frac{\langle S_{g}^{2}\rangle(a^{\prime})}{2}=$
$\displaystyle\frac{\langle
S_{g}^{2}\rangle(a)}{2}\left[1+\frac{2dg}{g}+(8K-4)dl\right]$
$\displaystyle+dl\frac{g^{2}}{2\pi^{3}}\int dr\mathrel{\mathop{\mathchar
58\relax}}(\nabla\phi(r))^{2}\mathrel{\mathop{\mathchar 58\relax}}.$ (49)

The last term $\langle S_{g}S_{\alpha}\rangle(a^{\prime})$ can be, yet again,
separated into a simple part and a part that requires the OPE 42 to give

$\displaystyle\langle S_{g}S_{\alpha}\rangle(a^{\prime})=$
$\displaystyle\langle
S_{g}S_{\alpha}\rangle(a)\Big{[}1+\frac{dg}{g}+\frac{d\alpha}{\alpha}$ (50)
$\displaystyle+(6K+s-4)dl\Big{]}+dl\frac{g}{\pi}\langle S_{\alpha}\rangle(a).$

Putting everything together, one arrives at

$\displaystyle\langle S_{g}\rangle(a^{\prime})+\langle
S_{g}\rangle(a^{\prime})-\frac{\langle
S_{g}^{2}\rangle}{2}(a^{\prime})-\langle S_{g}S_{\alpha}\rangle(a^{\prime})$
$\displaystyle=$ $\displaystyle\langle
S_{g}\rangle(a)\left[1+\frac{dg}{g}+(4K-2)dl-\frac{\alpha}{g}dl\right]$
$\displaystyle+\langle
S_{\alpha}\rangle(a)\left[1+\frac{d\alpha}{\alpha}+(2K+s-2)dl-\frac{g}{\pi}dl\right]$
$\displaystyle-\frac{\langle
S_{g}^{2}\rangle(a)}{2}\left[1+\frac{2dg}{g}+(8K-4)dl\right]$
$\displaystyle-\langle
S_{g}S_{\alpha}\rangle(a)\left[1+\frac{dg}{g}+\frac{d\alpha}{\alpha}+(6K+s-4)dl\right]$
$\displaystyle-dl\frac{\alpha(a)}{\pi^{2}u}\int
dxd\tau\langle\mathrel{\mathop{\mathchar
58\relax}}(\partial_{\tau}\phi(x,\tau))^{2}\mathrel{\mathop{\mathchar
58\relax}}\rangle$ $\displaystyle-dl\frac{g^{2}}{2\pi^{3}}\int
dr\mathrel{\mathop{\mathchar
58\relax}}(\nabla\phi(r))^{2}\mathrel{\mathop{\mathchar 58\relax}}.$ (51)

Upon imposing that the partition function $Z$ remains unchanged, the RG
equations for $\alpha$ and $g$ can be directly read off. Those for $K$ and $u$
are found by re-exponentiating the remaining quadratic terms. In the end, one
finds

$\displaystyle\frac{d}{dl}$ $\displaystyle\frac{u}{K}=\frac{g^{2}u}{\pi^{2}},$
(52) $\displaystyle\frac{d}{dl}$ $\displaystyle\frac{1}{uK}=\frac{2\alpha}{\pi
u}+\frac{g^{2}}{\pi^{2}u},$ (53) $\displaystyle\frac{d}{dl}$ $\displaystyle
g=(2-4K)g+\alpha,$ (54) $\displaystyle\frac{d}{dl}$
$\displaystyle\alpha=(2-s-2K)\alpha+\frac{g\alpha}{\pi}.$ (55)

##  Appendix C Gap equation from the variational method

We wish to find a solution to the self-consistent variational equation Eq. 26
using the following ansatz $G^{-1}_{\rm AFM}(q,\omega_{n})=\frac{1}{\pi
K}\left[uq^{2}+\frac{\omega_{n}^{2}}{u}+\frac{\Delta^{2}}{u}\right]$. Setting
$\omega_{n}=q=0$ in Eq. 26 yields

$\displaystyle\frac{\Delta^{2}}{\pi Ku}$
$\displaystyle=\frac{8g}{\pi^{2}a^{2}}\exp\left[-2K\int\frac{d\omega_{n}^{\prime}}{\sqrt{{\omega_{n}^{\prime}}^{2}+\Delta^{2}}}\right]$
(56) $\displaystyle+\frac{4\alpha
u^{-s}}{\pi^{2}a^{2-s}}\int_{\tau_{c}}^{\infty}\frac{d\tau}{\tau^{1+s}}\exp\left[-K\int
d\omega_{n}^{\prime}\frac{1+\cos\omega_{n}^{\prime}\tau}{\sqrt{{\omega_{n}^{\prime}}^{2}+\Delta^{2}}}\right],$

where we have done both integrals over $q$ in the exponentials. The remaining
integrals over $\omega_{n}^{\prime}$ are UV divergent and are regulated by
reintroducing the imaginary-time cutoff $\tau_{c}=a/u$ (thus $1/\tau_{c}$ in
Matsubara frequency space). In the limit of $\Delta\tau\ll 1$, one can then
find that
$\int_{0}^{1/\tau_{c}}\frac{d\omega_{n}^{\prime}}{\sqrt{{\omega_{n}^{\prime}}^{2}+\Delta^{2}}}=-\log(\Delta\tau_{c})$,
and
$\int_{0}^{1/\tau_{c}}d\omega_{n}^{\prime}\frac{1+\cos\omega_{n}^{\prime}\tau}{\sqrt{{\omega_{n}^{\prime}}^{2}+\Delta^{2}}}=-\log(\Delta^{2}\tau_{c}\tau)-\gamma_{E}$,
with $\gamma_{E}$ Euler’s gamma constant. Plugging this back into Eq. 56 gives

$\displaystyle\frac{\Delta^{2}}{\pi K}=$
$\displaystyle\frac{8g}{\pi^{2}\tau_{c}^{2-4K}}\Delta^{4K}$ (57)
$\displaystyle+\frac{4\alpha}{\pi^{2}\tau_{c}^{2-s-2K}}\int_{\tau_{c}}^{\frac{1}{\Delta}}\frac{d\tau}{\tau^{1+s-2K}}\Delta^{4K}e^{2K\gamma_{E}},$

which, after doing the integral over $\tau$ and rearranging the terms, leads
to

$\displaystyle\frac{(\Delta\tau_{c})^{2-4K}}{K}=$ $\displaystyle
g\frac{8}{\pi}+\alpha\frac{4e^{2K\gamma_{E}}}{\pi}\frac{(\Delta\tau_{c})^{s-2K}-1}{2K-s}.$
(58)

This is Eq. 28 in the text.

##  Appendix D Variational method vs RG

This section shows how the ansatz $G^{-1}_{\rm var}(q,\omega_{n})=\frac{1}{\pi
K_{r}}\left[u_{r}q^{2}+\frac{\omega_{n}^{2}}{u_{r}}+\frac{\Delta^{2}}{u_{r}}\right]$
fails at predicting the correct renormalized coefficients close to the
transition. For simplicity and to compare the result with the perturbative RG,
we focus on the LL phase where $\Delta_{r}=0$. Taking the second derivative
with respect to $q$ of Eq. 26 shows that $\frac{u_{r}}{K_{r}}=\frac{u}{K}$ so
we can express both Luttinger parameters in terms of $\eta$ such that
$K=K_{r}\sqrt{1+\eta}$, $u=u_{r}\sqrt{1+\eta}$. The ansatz thus becomes
$G^{-1}_{\rm var}(q,\omega_{n})=\frac{1}{\pi
K}\left[uq^{2}+\frac{\omega_{n}^{2}}{u}(1+\eta)\right]$ where $\eta$ is the
only parameter to solve for. Taking the second derivative of Eq. 26 with
respect to $\omega_{n}$ yields

$\displaystyle\eta=$ $\displaystyle\frac{\alpha
K}{\pi\tau_{c}^{2-s}}\int_{\tau_{c}}^{\infty}\frac{d\tau}{\tau^{-1+s}}$ (59)
$\displaystyle\times\exp\left[-\frac{K}{\pi}\int
dq^{\prime}d\omega_{n}^{\prime}\frac{1-\cos\omega_{n}^{\prime}\tau}{u{q^{\prime}}^{2}+\frac{{\omega_{n}^{\prime}}^{2}}{u}(1+\eta)}\right].$

We now perform the integral over $q$ and regulate the diverging integral over
$\omega_{n}^{\prime}$ by adding the imaginary-time UV cutoff $\tau_{c}$ to
obtain

$\displaystyle\eta=$ $\displaystyle\frac{\alpha
K}{\pi\tau_{c}^{2-s}}\int_{\tau_{c}}^{\infty}\frac{d\tau}{\tau^{-1+s}}$ (60)
$\displaystyle\times\exp\left[-2K_{r}\int_{0}^{1/\tau_{c}}d\omega_{n}^{\prime}\frac{1-\cos\omega_{n}^{\prime}\tau}{\omega_{n}^{\prime}}\right].$

In the limit of large $\tau/\tau_{c}$, one proves the following identity

$\displaystyle\int_{0}^{1/\tau_{c}}d\omega_{n}^{\prime}\frac{1-\cos\omega_{n}^{\prime}\tau}{\omega_{n}^{\prime}}=\log(\tau/\tau_{c})+\gamma_{E}+\mathcal{O}(\tau_{c}/\tau),$
(61)

which implies that

$\displaystyle\eta=$ $\displaystyle\frac{\alpha
Ke^{-2K_{r}\gamma_{E}}}{\pi\tau_{c}^{2-s-2K_{r}}}\int_{\tau_{c}}^{\infty}\frac{d\tau}{\tau^{-1+s+2K_{r}}}.$
(62)

Computing this integral and using $K=K_{r}\sqrt{1+\eta}$ leads to

$\displaystyle\frac{\eta}{\sqrt{1+\eta}}=$ $\displaystyle\frac{\alpha
K_{r}e^{-2K_{r}\gamma_{E}}}{\pi(2K_{r}+s-2)}.$ (63)

Now comes the problem. From the argument made in section V, it is clear that
the transition is given in terms of $K_{r}$ by $(K_{r})_{c}=\max(1-s/2,1/2)$.
This means that for $s<1$, Eq. 63 predicts a renormalization of $K$ and $u$
that diverges towards the transition as

$\displaystyle\sqrt{\eta}\simeq$ $\displaystyle\frac{\alpha
K_{r}e^{-2K_{r}\gamma_{E}}}{\pi(2K_{r}+s-2)}.$ (64)

which is of course highly unphysical. This result can however be recovered
from the RG using a very crude approximation: linearizing the RG flow about
$\alpha=0$ (see fig. 6). To show how this recovers Eq. 64, let us start by
recalling the important RG equations for $s<1$ (see subsection IV.3)

Figure 6: RG flow in the LL phase of the couplings $\alpha$ and $K$ according
to the actual RG equations (solid gold lines) and the linearized equations
(dashed red lines). This plot is for $s=0.5$.

$\displaystyle\frac{d}{dl}$ $\displaystyle K=-\frac{\alpha K^{2}}{\pi},$ (65)
$\displaystyle\frac{d}{dl}$ $\displaystyle\alpha=(2-s-2K)\alpha.$ (66)

In the LL phase, $\alpha$ flows to 0 and $K$ to $K_{r}$. This means that
$K_{r}$ is given by
$K_{r}=K+\int_{K}^{K_{r}}dK^{\prime}=K+\int_{\alpha}^{0}\frac{dK}{d\alpha}(\alpha^{\prime})d\alpha^{\prime}$.
For the linearized RG flow shown in fig. 6,
$\frac{dK}{d\alpha}(\alpha^{\prime})=\frac{dK}{d\alpha}(\alpha^{\prime}=0)=\frac{K_{r}^{2}}{\pi(2K_{r}+s-2)}$.
Plugging this back into Eq. 66 and using $K=K_{r}\sqrt{1+\eta}$ yields

$\displaystyle\sqrt{1+\eta}-1=$ $\displaystyle\frac{\alpha
K_{r}}{\pi(2K_{r}+s-2)}.$ (67)

Figure 7: Numerical evaluation of observables for the commensurate dissipative
spin chain. The values of the parameters are taken to be $s=1$ ($K_{c}=0.5$),
$K=0.65$, $g=0.5$, $u=1$, and $\beta$ is scaled as $L$. For the left and
middle column, red, green, and blue colors denote $L=128$, $L=256$ and $L=320$
respectively. All the quantities have been averaged over 10000-20000
configurations and the value of time-step $dt$ is $0.5$. (top) and (bottom)
rows correspond to $\alpha=0.5$ (LL) and $\alpha=6$ (AFM), respectively.
(left) Susceptibility $\chi$ as a function $q$. For $\alpha=0.5$, it stays
finite and constant whereas for $\alpha=6$, it vanishes for small $q$.
(middle) Behavior of $C(\omega_{n})$ as a function of $\omega_{n}$. For
$\alpha=0.5$, $\omega_{n}C(\omega_{n})$ saturates to a constant
$K_{r}/2=0.307$; for $\alpha=6$, $C(\omega_{n})$ becomes a constant
$K_{r}/(2\Delta_{r})=0.07$. (right) Behavior of order parameter
$\langle\cos\left[2\left(\phi-\phi_{\rm CoM}\right)\right]\rangle$ as a
function of system size $L$. For $\alpha=0.5$, the order parameter decays as
$L^{-K_{r}}$ with $K_{r}=0.59$, whereas for $\alpha=6$, it increases and
saturates to a constant as $c_{1}-c_{2}/L$ with $c_{1}=0.742$ and
$c_{2}=0.0296$.

Near the transition ($K_{r}\to(1-s/2)^{+}$), this equation becomes exactly Eq.
64 up to the numerical factor $e^{-2K_{r}\gamma_{E}}$ which is probably an
artefact of the approximation made in Eq. 61. The main takeaway message is
that the variational method captures the RG flow linearized about $\alpha=0$,
and thus breaks down close to the transition.

Figure 8: Numerical evaluation of observables for the commensurate dissipative
spin chain. The values of the parameters are taken to be $s=1.5$
($K_{c}=0.5$), $K=0.75$, $g=1$, $u=1$, and $\beta$ is scaled as $L$. For the
left and middle column, red, green, and blue colors denote $L=128$, $L=256$
and $L=320$ respectively, and $\beta$ is scaled as $L$. All the quantities
have been averaged over 10000-20000 configurations and the value of time-step
$dt$ is $0.5$. (top) and (bottom) rows correspond to $\alpha=0.65$ (LL) and
$\alpha=6$ (AFM), respectively. (left) Susceptibility $\chi$ as a function
$q$. For $\alpha=0.65$, it stays finite and constant whereas for $\alpha=6$,
it vanishes for small $q$. (middle) Behavior $C(\omega_{n})$ as a function of
$\omega_{n}$. For $\alpha=0.65$, $\omega_{n}C(\omega_{n})$ saturates to a
constant $K_{r}/2=0.345$; for $\alpha=6$, $C(\omega_{n})$ becomes a constant
$K_{r}/(2\Delta_{r})=0.031$. (right) Behavior of order parameter
$\langle\cos\left[2\left(\phi-\phi_{\rm CoM}\right)\right]\rangle$ as a
function of system size $L$. For $\alpha=0.65$, the order parameter decays as
$L^{-K_{r}}$ with $K_{r}=0.637$, whereas for $\alpha=5$, it increases and
saturates to a constant as $c_{1}-(c_{2}/L)$ with $c_{1}=0.872$ and $c_{2}=0$.

##  Appendix E A simple argument

In this section we show that the previous arguments allow to correctly capture
the phase transition if the renormalisation of $K$ is taken into account in a
simpler manner, namely without resorting to a self consistent computation.
This approximation should correspond to assuming that all the way down to
$K_{r}=K_{c}$ the computation of $K_{r}$ is perturbative in $\alpha$. The
computation amounts to replacing, in Eq. 62, the value of $K_{r}$ with $K$ as

$\eta=\frac{\alpha KC}{2\pi(K-K_{c})},$ (68)

with $C=e^{-2K\gamma_{E}}$, and then using $\eta$’s definition

$\eta=\left(\frac{K}{K_{r}}\right)^{2}-1,$ (69)

which, for $K_{r}=K_{c}$, recovers

$\left(K-K_{c}\right)^{2}\propto\alpha,$ (70)

which is the parabolic shape of the transition predicted by RG.

##  Appendix F Derivation of the order parameter

Figure 9: $K_{r}/(2\Delta_{r})-C(\omega_{n})$ for the dissipative phase
($s=1.5,K=0.75,g=1,\alpha=8$). This quantity fits well with $9\times
10^{-5}\omega_{n}^{2}$ for small $\omega_{n}$, which is the sub-leading
$\omega_{n}$ dependence of the propagator in the AFM.

In this appendix, we show the calculation of the order parameter
$\langle\cos\left(2\left(\phi-\phi_{\rm CoM}\right)\right)\rangle$ using our
Gaussian variational ansatz. The behavior of this quantity in the LL has
already been calculated in [13] (appendix A shows that
$\langle\cos\left(2\left(\phi-\phi_{\rm CoM}\right)\right)\rangle\propto
L^{-K_{r}}$) so we concentrate on the AFM phase defined by the propagator
$G_{\rm AFM}^{-1}(q,\omega_{n})=\frac{1}{\pi
K}(uq^{2}+\frac{\omega_{n}^{2}}{u}+\frac{\Delta_{r}^{2}}{u_{r}})$. Since we
work within our Gaussian variational theory, one writes $\langle\cos
2\left(\phi-\phi_{\rm CoM}\right)\rangle=\exp(-2\langle(\phi-\phi_{\rm
CoM})^{2}\rangle)$ where $\langle(\phi-\phi_{\rm CoM})^{2}\rangle$ can be
broken up as

$\displaystyle\langle(\phi-\phi_{\rm CoM})^{2}\rangle=\frac{1}{\beta
L}\left[\sum\limits_{q\neq 0}G_{\rm AFM}(q,0)\right.$ $\displaystyle+$
$\displaystyle\sum\limits_{\omega_{n}\neq 0}G_{\rm
AFM}(0,\omega_{n})+\left.\sum\limits_{\omega_{n}\neq 0,q\neq 0}G_{\rm
AFM}(q,\omega_{n})\right].$ (71)

The first two terms can be computed for finite $\beta$ and $L$ using Matsubara
sum techniques. This leads to

$\displaystyle\sum_{q\neq 0}G_{\text{AFM}}(q,0)=$ $\displaystyle\frac{\pi
K_{r}L}{2\Delta_{r}}\coth\left(\frac{\Delta_{r}L}{2u_{r}}\right)-\frac{\pi
K_{r}u_{r}}{\Delta_{r}^{2}},$ (72)

and

$\displaystyle\sum_{\omega_{n}\neq 0}G_{\text{AFM}}(0,\omega_{n})=$
$\displaystyle\frac{\pi
K_{r}u_{r}\beta}{2\Delta_{r}}\coth\left(\frac{\Delta_{r}\beta}{2}\right)-\frac{\pi
K_{r}u_{r}}{\Delta_{r}^{2}}.$ (73)

The third term is then replaced by its continuous limit
($L\to\infty,\beta\to\infty$) as the finite size effects are taken care of in
the first two sums. Thus

$\displaystyle\frac{1}{\beta L}$ $\displaystyle\sum_{\omega_{n}\neq 0,q\neq
0}G_{\rm AFM}(q,\omega_{n})$ $\displaystyle=$
$\displaystyle\frac{K_{r}}{4\pi}\underset{\sqrt{(u_{r}q)^{2}+\omega_{n}^{2}}<\frac{1}{\tau_{c}}}{\int
d\omega_{n}dq}\frac{1}{u_{r}q^{2}+\frac{\omega_{n}^{2}}{u_{r}}+\frac{\Delta_{r}^{2}}{u_{r}}}$
$\displaystyle=$
$\displaystyle-\frac{K_{r}}{4}\log\left(\frac{(\tau_{c}\Delta)^{2}}{1+(\tau_{c}\Delta)^{2}}\right),$
(74)

where we have introduced the UV cutoff $\tau_{c}$. In our simulations, we take
$L=\beta$ and then $L\to\infty$. Taking this same scaling, a large $L$
expansion of the order parameter reads

$\displaystyle\langle\cos(2(\phi-\phi_{\rm CoM}))\rangle\simeq$ $\displaystyle
c_{1}-\frac{c_{2}}{L}+\frac{c_{3}}{L^{2}},$ (75)

where
$c_{1}=\left(\frac{\Delta\tau_{c}}{\sqrt{1+(\Delta\tau_{c})^{2}}}\right)^{K_{r}}$,
$c_{2}=c_{1}\frac{\pi K_{r}}{\Delta_{r}}\left(1+u_{r}\right)$, and
$c_{3}=c_{1}\frac{\pi K_{r}}{\Delta_{r}^{2}}\left[4u_{r}+\pi
K_{r}\left(1+u_{r}\right)^{2}\right]$ respectively. We observe that the order
parameter monotonically increases to a constant with a leading finite size
dependence of $1/L$. This is a signature of the gapped phase, as a fractional
phase would have finite size scaling $\propto L^{\frac{s}{2}-1}$, and a
Luttinger liquid would have a $\propto L^{-K_{r}}$ behaviour.

##  Appendix G Additional numerical results

In this appendix, we provide additional numerical results that fortify some of
the analytical claims made in the main article. Firstly, we calculate the
correlation functions discussed in section VI for ohmic and superohmic baths.
We then provide numerical evidence of the absence of the fractional excitation
$|\omega_{n}|^{s}$ term in the gapped phase. The results are provided in fig.
7 and fig. 8, where numerical results are given for ohmic ($s=1$) and
superohmic ($s=1.5$) baths similarly to the main article.

We also provide numerical evidence of the absence of a fractional excitation
term $|\omega_{n}|^{s}$ in fig. 9. One can compute $C(\omega_{n})$ using the
AFM variational propagator $G_{\rm var}=\pi
K_{r}\left[uq^{2}+\frac{\omega_{n}^{2}}{u}+\frac{\Delta_{r}^{2}}{u_{r}}\right]^{-1}$
and see that
$C(\omega_{n})=\frac{K_{r}}{2}\frac{1}{\sqrt{\omega_{n}^{2}+\Delta_{r}^{2}}}\simeq\frac{K_{r}}{2\Delta_{r}}\left[1-\frac{1}{2}\left(\frac{\omega_{n}^{2}}{\Delta_{r}}\right)^{2}\right]$.
Thus, one can subtract $C(\omega_{n})$ from the constant $K_{r}/(2\Delta_{r})$
and plot the quantity as a function of $\omega_{n}$ in log-log scale to check
the sub-leading $\omega_{n}$-dependence. This is done in fig. 9, and we find
that the exponent is approximately $2$ as expected.

## References

  * Anderson [1958] P. W. Anderson, Phys. Rev. 109, 1492 (1958).
  * Maimbourg  _et al._ [2021] T. Maimbourg, D. M. Basko, M. Holzmann, and A. Rosso, Phys. Rev. Lett. 126, 120603 (2021).
  * Caldeira and Leggett [1983] A. Caldeira and A. Leggett, Annals of Physics 149, 374 (1983).
  * Bray and Moore [1982] A. J. Bray and M. A. Moore, Phys. Rev. Lett. 49, 1545 (1982).
  * Leggett  _et al._ [1987] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59, 1 (1987).
  * Fisher and Zwerger [1985] M. P. A. Fisher and W. Zwerger, Phys. Rev. B 32, 6190 (1985).
  * Schmid [1983] A. Schmid, Phys. Rev. Lett. 51, 1506 (1983).
  * Nandkishore and Huse [2015] R. Nandkishore and D. A. Huse, Annual Review of Condensed Matter Physics 6, 15 (2015).
  * Altshuler  _et al._ [1997] B. L. Altshuler, Y. Gefen, A. Kamenev, and L. S. Levitov, Phys. Rev. Lett. 78, 2803 (1997).
  * Abanin  _et al._ [2019] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev. Mod. Phys. 91, 021001 (2019).
  * De Filippis  _et al._ [2020] G. De Filippis, A. de Candia, L. M. Cangemi, M. Sassetti, R. Fazio, and V. Cataudella, Physical Review B 101, 10.1103/physrevb.101.180408 (2020).
  * Majumdar  _et al._ [2023a] S. Majumdar, L. Foini, T. Giamarchi, and A. Rosso, Phys. Rev. B 107, 165113 (2023a).
  * Majumdar  _et al._ [2023b] S. Majumdar, L. Foini, T. Giamarchi, and A. Rosso, Phys. Rev. B 108, 205138 (2023b).
  * Cazalilla  _et al._ [2006] M. A. Cazalilla, F. Sols, and F. Guinea, Physical Review Letters 97, 10.1103/physrevlett.97.076401 (2006).
  * Cai  _et al._ [2014] Z. Cai, U. Schollwöck, and L. Pollet, Physical Review Letters 113, 10.1103/physrevlett.113.260403 (2014).
  * Malatsetxebarria  _et al._ [2013] E. Malatsetxebarria, Z. Cai, U. Schollwöck, and M. A. Cazalilla, Physical Review A 88, 10.1103/physreva.88.063630 (2013).
  * Giamarchi and Schulz [1987] T. Giamarchi and H. J. Schulz, Europhysics Letters 3, 1287 (1987).
  * Giamarchi and Schulz [1988] T. Giamarchi and H. J. Schulz, Phys. Rev. B 37, 325 (1988).
  * Giamarchi [2004] T. Giamarchi,  _Quantum physics in one dimension_, International series of monographs on physics (Clarendon Press, Oxford, 2004).
  * von Delft and Schoeller [1998] J. von Delft and H. Schoeller, Annalen der Physik 510, 225 (1998).
  * Giachetti  _et al._ [2021] G. Giachetti, N. Defenu, S. Ruffo, and A. Trombettoni, Physical Review Letters 127, 10.1103/physrevlett.127.156801 (2021).
  * Giachetti  _et al._ [2022] G. Giachetti, A. Trombettoni, S. Ruffo, and N. Defenu, Physical Review B 106, 10.1103/physrevb.106.014106 (2022).
  * Ribeiro  _et al._ [2023] A. L. S. Ribeiro, P. McClarty, P. Ribeiro, and M. Weber, Dissipation-induced long-range order in the one-dimensional bose-hubbard model (2023), arXiv:2311.07683 [cond-mat.str-el] .
  * Weber [2023] M. Weber, Quantum spin chains with bond dissipation (2023), arXiv:2310.11525 [cond-mat.str-el] .
  * Note [1] Here and in the following, we call BKT transition a transition whose RG equations are of the BKT type.
  * Cardy [1996] J. Cardy,  _Scaling and Renormalization in Statistical Physics_, Cambridge Lecture Notes in Physics (Cambridge University Press, 1996).
  * Dupuis [2023] N. Dupuis,  _Field Theory of Condensed Matter and Ultracold Gases_ (WORLD SCIENTIFIC (EUROPE), 2023).
  * Note [2] Pushing the RG equations to higher orders shows that $\frac{d}{dl}\frac{u}{K}=\mathcal{O}(\alpha^{2})$. This is to be contrasted with the incommensurate spin chain [12, 13] for which $\frac{d}{dl}\frac{u}{K}=0$ is non-perturbative. This is due to the fact that the system-bath interaction in the incommensurate spin chain action is invariant under the statistical tilt symmetry $\phi(x,\tau)\to\phi(x,\tau)+\mu x$, where $\mu$ is a constant.
  * Maghrebi  _et al._ [2017] M. F. Maghrebi, Z.-X. Gong, and A. V. Gorshkov, Physical Review Letters 119, 10.1103/physrevlett.119.023001 (2017).
  * Feynman [1988] R. Feynman,  _Statistical Mechanics: A Set of Lectures_, Frontiers in physics (Sarat Book Distributors, 1988).
  * Rajaraman [1982] R. Rajaraman,  _Solitons and Instantons: An Introduction to Solitons and Instantons in Quantum Field Theory_, North-Holland personal library (North-Holland Publishing Company, 1982).
  * Maki [1978] K. Maki, Phys. Rev. B 18, 1641 (1978).
  * HIDA and NAKAYA [1986] K. HIDA and S. NAKAYA, in  _Proceedings of the Yamada Conference XV on Physics and Chemistry of Quasi One-Dimensional Conductors_, edited by S. Tanaka and K. Uchinokura (Elsevier, 1986) pp. 140–142.
  * Huse  _et al._ [2015] D. A. Huse, R. Nandkishore, F. Pietracaprina, V. Ros, and A. Scardicchio, Phys. Rev. B 92, 014203 (2015).
  * Minnhagen [1987] P. Minnhagen, Rev. Mod. Phys. 59, 1001 (1987).


# The Concept of Forward-Forward Learning Applied to a Multi Output Perceptron

K. F. Karlsson Correspondence to<EMAIL_ADDRESS>Saab
Dynamics, SE-581 88 Linköping, Sweden

(06 April 2023)

###### Abstract

The concept of a recently proposed Forward-Forward learning algorithm for
fully connected artificial neural networks is applied to a single multi output
perceptron for classification. The parameters of the system are trained with
respect to increased (decreased) ”goodness” for correctly (incorrectly)
labelled input samples. Basic numerical tests demonstrate that the trained
perceptron effectively deals with data sets that have non-linear decision
boundaries. Moreover, the overall performance is comparable to more complex
neural networks with hidden layers. The benefit of the approach presented here
is that it only involves a single matrix multiplication.

##

The perceptron is a type of artificial neural network that was introduced in
the late 1950s by psychologist Frank Rosenblatt as a model of a biological
neuron [1]. It is a single-layer feedforward neural network that learns from
examples to classify input patterns by adjusting the weights of its
connections. The perceptron has been widely used in pattern recognition and
machine learning applications due to its simplicity and efficiency.

Despite its early success, the perceptron has a major limitation: it can only
learn to distinguish linearly separable patterns [2]. This means that it can
only classify patterns that can be separated by a hyperplane in the input
space. Patterns that are not linearly separable, such as the Boolean exclusive
OR (XOR) function, cannot be learned by a single perceptron.

To overcome this limitation, multiple-layer neural networks were developed,
such as the multilayer perceptron (MLP), which consists of several layers of
interconnected neurons where the input to an interior hidden layer is given by
the output of the preceding layer. These networks can learn non-linear
decisions boundaries that separate the patterns, and it has been shown that
nearly any function, to arbitrary accuracy, can be approximated by such a
network with one hidden layer of sufficiently many neurons [3].

Supervised training of neural networks is today typically accomplished by
backpropagation, a learning algorithm that adjusts the weights of the
connections based on the error between the actual output and the desired
output. Backpropagation works by propagating the error backwards through the
network, layer by layer, allowing all its weights to be adjusted in a way that
minimises the error [4].

For a typical multiclass classification problem, the input dimension of the
perceptron is given by the dimension of the input signal, and the dimension of
the output vector is equal to the number of classes. For a given input, the
output value for each dimension can be interpretable by the softmax function
as the probability for the input signal that it belongs to the corresponding
class.

Recently, G. Hinton proposed the Forward-Forward (FF) algorithm that instead
of minimising the final error at the network output optimizes the ”goodness”
of each involved network layer individually [5]. Thus, this algorithm does not
rely on backpropagation of information to preceding layers during learning.
The goodness of a layer is a measure of the length of the output vector, and
can be defined as the summed square of its elements. The aim is that good or
bad input data should yield high or low goodness, respectively (referred to by
Hinton as real and negative data). In this context, the procedure of
supervised learning of a classifier is to provide the input not only with the
training signal, but also with the the actual class label. If more than two
classes are involved, one hot encoding is employed. As a consequence, the
input vector is extended by one extra dimension per class. Input data with
correct labels correspond to good data, and input data with incorrect labels
correspond to bad data.

In the above-mentioned publication, a network with four hidden layers of 2000
neurons each were reported to yield a performance comparable to state-of-art
networks conventionally trained by backpropagation [5]. Only the activity
vectors of the last three hidden layers are combined in order to extract a
collective goodness value, whereas the goodness of the first layer is omitted.
By excluding the first layer, the overall performance reported to increase.

In the present work, a single perceptron trained and used within the FF
framework is studied. The training is performed with the Adam variant of
stochastic gradient descent [6], to minimise the logistic loss applied to the
goodness value minus a threshold. First, it will be shown that a single FF
perceptron can be trained to mimic the XOR function. In this case, at least
three inputs are required for training, the two XOR inputs and the class
label. After training, the perceptron is expected to yield high goodness only
if the two XOR input values correspond to the correct label, also given as
input. Otherwise, low goodness is expected. Thus, this perceptron constitutes
a binary classifier with the two categories of correctly and incorrectly
labelled inputs. Since these two categories are not linearly separable in the
three-dimensional input space, the expected goodness cannot be achieved by a
single output. By using multiple outputs, however, the goodness value is
determined by the collective properties of all outputs and this enables non-
linear decision boundaries. As a consequence, a single FF trained perceptron,
with more than one output, can perform XOR classification (see Fig. 1 for some
examples with different number of outputs and activation functions). Thus,
adding outputs to this perceptron is analogous to adding neurons to a hidden
layer in a conventional MLP.

Figure 1: Classifications maps of a perceptron with two (a and b) or four (c)
outputs trained to mimic the XOR function. Boolean true and false are
represented by values near -1 or near +1 at the two inputs $x_{1}$ and
$x_{2}$, respectively. The two classes are labelled by 0 for equal inputs, and
1 for different inputs. Given $x_{1}$ and $x_{2}$, the classification is based
on the label that yields highest goodness, a value defined as the sum of the
squared output vector elements $\mathbf{y}$. Above each map is the
corresponding mathematical expression of the perceptron, including values for
the trained parameters.

A slightly more complex two-dimensional example is based on classifying noisy
input points $(x_{1},x_{2})$ scattered along two blue and red spiralling
curves according to Fig. 2a. Points of different colour belong to different
classes. Besides $x_{1}$ and $x_{2}$, the class label should also be provided
at the input. It is here given at a single input, label 0 for blue and label 1
for red, elevating the input dimension to three. With 32 leaky $ReLU$
activated outputs, the perceptron has 128 trainable parameters (including 32
biases). It classifies more than 99 % of random test points correctly after
training (see Fig. 2b).

Figure 2:  Training data scattered around spiralling curves (a) of blue and
red points belonging to two classes labelled 0 and 1, respectively, and (b)
the resulting classification map of a trained perceptron with 3 inputs and 32
outputs. More than 99% of random tests points are correctly classified by the
trained perceptron.

The previous examples clearly demonstrate that a single multi-output
perceptron, trained and used within the framework of FF, is able to handle
data sets with complex decision boundaries. Moreover, the number of parameters
can be made arbitrary large by just increasing the number of outputs, and the
performance is expected to improve with more parameters. As a final basic
example, with input and output spaces of significantly higher dimensions, the
perceptron with $ReLU$ activated outputs was trained to classify images of
handwritten digits using the famous MNIST dataset [7]. 60000 images were used
for training and 10000 different images were used for evaluating the trained
perceptron. The training set was split into mini-batches of size 10, and the
training was performed such that 10 training image first was provided with
correct labels (one-hot encoded), as good data, and subsequently the same
images were provided with a random incorrect label as bad data. Gradients were
computed and weights were updated after each mini-batch. See the code for more
details [8].

Results obtained for some selected numbers of outputs are summarised in Table
1, with a test error below 1.7 % for 2000 or more outputs. Augmenting the
training data by randomly jittering the images up to two pixels in any
direction, reduces the test error below 1.0% for the perceptron with 8000
outputs.

Outputs | Data Augmentation | Test Error $\lesssim$ (%) | Epochs
---|---|---|---
125 | none | 2.6 | 80
500 | none | 1.9 | 80
2000 | none | 1.7 | 150
8000 | none | 1.6 | 150
125 | jittering | 2.2 | 180
500 | jittering | 1.4 | 260
2000 | jittering | 1.1 | 340
8000 | jittering | 1.0 | 340
Table 1: Results obtained on classification of the MNIST validation data after
FF training of a single perceptron with various number of outputs. After the
listed number of epochs, no significant improvement of the test error could be
observed. A learning rate of 0.0003 was used, the goodness was computed as the
mean of the squared output elements and threshold in the logistic function was
set to 10 [8].

In many respects, the single FF perceptron studied here behaves as a
conventional two-layer perceptron, i.e., with one hidden layer. It is
considerably slower to train, however, and it seems to adapt less accurately
to the training data than the conventional MLP. This property may make the FF
perceptron less prone to overfitting.

In the future, it is expected that analogue neural networks will play an
important role, mainly due to their predicted speed and low energy
consumption. The single matrix multiplication involved in the evaluation of
the perceptron outputs, combined with the absence of back propagation of
errors between layers, makes the hardware architecture for analogue matrix
computations much simpler. Indeed, single trainable perceptrons have already
been realized in integrated memristor chips [9], although still with very few
inputs and outputs. Thus, high performing single perceptrons, as the one
studied here, may find use for fast and energy efficient real time signal
classification by future analogue hardware.

Artificial neural networks are often regarded as black-box systems, which can
be problematic in critical domains such as medical or military applications
due to their lack of explainability. However, a perceptron that uses a single
matrix for input-output mapping provides a transparent link between inputs and
outputs, which promotes interpretability and explainability of the network.
This transparency makes it easier to both visualize and understand the factors
that influence outputs of the network.

In conclusion, it has been demonstrated that the Forward-Forward learning
procedure enables a single perceptron to be trained to mimic the XOR function
and to handle more complex and non-linear decision boundaries. For basic tests
made on the famous MNIST data set of handwritten digits, the performance is
comparable to conventional neural networks with hidden layers, and a
classification test error smaller than 1% was achieved. The simplicity of a
single matrix multiplication promotes explainablity, and combined with the
absence of backpropagation it opens up for simpler analogue computational
hardware.

## Acknowledgements

The author is thankful to C. Larsson for bringing the recently proposed
Forward-Forward algorithm to his attention, and to A. Gällström for
interesting discussions on this topic.

## References

  * [1] Rosenblatt, F. The perceptron: A probabilistic model for information storage and organization in the brain. Psychological Review. 65, 386-408 (1958), http://doi.org/10.1037/h0042519
  * [2] Minsky, M. & Papert, S. Perceptrons. (MIT Press,1969)
  * [3] Cybenko, G. Approximation by superpositions of a sigmoidal function. Mathematics Of Control, Signals, And Systems (MCSS). 2, 303-314 (1989,12), http://doi.org/10.1007/BF02551274
  * [4] Rumelhart, D., Hinton, G. & Williams, R. Learning representations by back-propagating errors. Nature. 323, 533-536 (1986), http://doi.org/10.1038/323533a0
  * [5] Hinton, G. The Forward-Forward Algorithm Some Preliminary Investigations. arXiv:2212.13345 (2022), http://doi.org/10.48550/arXiv.2212.13345
  * [6] Kingma, D. & Ba, J. Adam: A Method for Stochastic Optimization. 3rd International Conference On Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings. (2015)
  * [7] Deng, L. The MNIST database of handwritten digit images for machine learning research. IEEE Signal Processing Magazine. 29, 141-142 (2012)
  * [8] K. F. Karlsson, Forward-Forward Learning of Perceptron. GitHub Repository. (2023), http://github.com/KFKarlsson/FFPerceptron
  * [9] Cai, F., Correll, J., Lee, S., Lim, Y., Bothra, V., Zhang, Z., Flynn, M. & Lu, W. A fully integrated reprogrammable memristor–CMOS system for efficient multiply–accumulate operations. Nature Electronics. 2 pp. 1 (2019,7), http://doi.org/10.1038/s41928-019-0270-x