aainstitutetext: Université Paris-Saclay, CNRS, LPTMS, 91405, Orsay,
Francebbinstitutetext: Univ Lyon, ENS de Lyon, CNRS, Laboratoire de Physique,
F-69342 Lyon, France

# Entanglement entropy of two disjoint intervals
and spin structures in interacting chains
in and out of equilibrium

Vanja Marić b Saverio Bocini a Maurizio Fagotti

###### Abstract

We take the paradigm of interacting spin chains, the Heisenberg
spin-$\frac{1}{2}$ XXZ model, as a reference system and consider interacting
models that are related to it by Jordan-Wigner transformations and
restrictions to sub-chains. An example is the fermionic analogue of the
gapless XXZ Hamiltonian, which, in a continuum scaling limit, is described by
the massless Thirring model. We work out the Rényi-$\alpha$ entropies of
disjoint blocks in the ground state and extract the universal scaling
functions describing the Rényi-$\alpha$ tripartite information in the limit of
infinite lengths. We consider also the von Neumann entropy, but only in the
limit of large distance. We show how to use the entropies of spin blocks to
unveil the spin structures of the underlying massless Thirring model. Finally,
we speculate about the tripartite information after global quenches and
conjecture its asymptotic behaviour in the limit of infinite time and small
quench. The resulting conjecture for the “residual tripartite information”,
which corresponds to the limit in which the intervals’ lengths are infinitely
larger than their (large) distance, supports the claim of universality
recently made studying noninteracting spin chains. Our mild assumptions imply
that the residual tripartite information after a small quench of the
anisotropy in the gapless phase of XXZ is equal to $-\log 2$.

##  1 Introduction

The entanglement entropies of subsystems in critical chains at low temperature
have become standard tools to extract universal properties of the underlying
conformal field theories Amico2008 ; Calabrese2009Entanglement . To the best
of our knowledge, the first evidence of it dates back to 1994, when Ref.
Holzhey1994Geometric  has shown that the leading order of the asymptotic
expansion of the von Neumann entropy of an interval in the limit of large
length is completely characterised by the central charge $c$

$S[A]=-\mathrm{tr}[\rho_{A}\log\rho_{A}]\sim\frac{c}{3}\log|A|\,.$ (1.1)

Here $\rho_{A}$ is the reduced density matrix of the interval $A$ and $|A|$
denotes $A$’s length. About ten years later that result was generalised
Calabrese2004Entanglement  to other settings, e.g. finite temperature and
finite systems. A similar formula was also obtained for the Rényi entropies

$S_{\alpha}[A]=\frac{1}{1-\alpha}\log\mathrm{tr}[\rho_{A}^{\alpha}]\sim\frac{c}{6}\left(1+\frac{1}{\alpha}\right)\log|A|\,.$
(1.2)

The interest in the subject exploded. On the one hand, it was found that some
corrections to the scaling are universal and determined by the conformal
dimensions of the operators of the theory Calabrese2010Parity ; Cardy2010 . On
the other hand, it was realised that changing the topology of the subsystem
uncovers additional universal properties Caraglio2008Entanglement ;
Furukawa2009Mutual . A huge effort was then directed to the study of the
entanglement entropies of disjoint blocks Casini2009Remarks ; Rajabpour2012 ;
Calabrese2009Entanglement1 ; Calabrese2011Entanglement ; Coser2014OnRenyi ;
Ruggiero2018 ; Alba2010Entanglement ; Blanco2011 ; Alba2011Entanglement ;
Fagotti2010disjoint ; Fagotti2012New ; Fries2019 ; Balasubramanian2011 ;
Grava2021 ; Ares2021Crossing . The problem was traced back to the asymptotic
behaviour of partition functions on Riemann surfaces with arbitrarily high
genus, and few explicit solutions have been obtained. We mention the Ising
model Calabrese2011Entanglement  and the free boson compactified on a circle
Calabrese2009Entanglement1  or on an orbifold Alba2011Entanglement .
Restricting to two disjoint blocks, a key role is played by a linear
combination of the entropies of subsystems, built up with the two blocks and
the region in between, which is also known as tripartite information

$I_{3}^{(\alpha)}(A,B,C)=S_{\alpha}[A]+S_{\alpha}[B]+S_{\alpha}[C]-S_{\alpha}[A\cup
B]-S_{\alpha}[A\cup C]-S_{\alpha}[B\cup C]+S_{\alpha}[A\cup B\cup C]\,.$ (1.3)

By construction, that combination simplifies both contributions coming from
the boundaries of the subsystems and those extensive with subsystems’ lengths,
and it turns out that it is even free from the log divergences (in the
subsystems’ lengths) that characterise the entanglement entropies at
criticality Casini2009Remarks ; Hayden2013Holographic . In the ground state of
a conformal field theory, $I_{3}^{(\alpha)}(A,B,C)$ turns out to be a
nontrivial function of the cross ratio

$x=\frac{|A||C|}{(|A|+|B|)(|B|+|C|)}$ (1.4)

and exhibits the symmetry $x\leftrightarrow 1-x$.

An important aspect that distinguishes the entropy of disjoints blocks from
that of an interval is its sensitivity to duality mappings, in that such
transformations generally preserve the locality of only a subset of operators
(including the Hamiltonian density), breaking it instead into semilocality for
the rest of them. This affects the meaning itself of subsystem, especially
when it is not simply connected. A consequence of it on the behaviour of the
entanglement entropies was firstly observed in the ground state of XY spin
chains Igloi2010 ; Fagotti2010disjoint . At the critical point the tripartite
information has a nontrivial scaling limit despite vanishing in the
corresponding free-fermion model that emerges after a Jordan-Wigner
transformation Casini2005 ; Casini2009reduced_density ; Casini2009Entanglement
.

Since the first step towards an analytical investigation into the entanglement
entropies of spin blocks is often a mapping to quasiparticles, it is important
to understand how such mappings affect the entropies and how to ultimately fix
the lack of equivalence. In that respect, Ref. Fagotti2010disjoint  proposed
an expansion of the tripartite information based on the decomposition of the
reduced density matrix into four pseudo reduced density matrices (pseudo-RMDs)
that unfold the spin-flip symmetry. On the other hand, the critical spin-chain
models studied in Ref. Fagotti2010disjoint  admit a quantum field theory
description of the low-energy and large-distance limit (the Ising model is
described by a massless Majorana field theory, while the XX model by a Dirac
one); as mentioned above, the Rényi entropy $S_{\alpha}$ of two disjoint
intervals can then be written in terms of partition functions on particular
Riemann surfaces. Ref. Coser2016Spin  examined what the lattice procedure of
Ref. Fagotti2010disjoint  translates into in the quantum field theory
description and established a one-to-one correspondence between the partition
functions with given spin structures and the moments of the pseudo-RDMs
introduced in Ref. Fagotti2010disjoint . It was argued, in particular, that
the entropy of disjoint blocks of spins is obtained by summing over all spin
structures of the fermion model on the underlying Riemann surface, i.e., over
all combinations of antiperiodic, also called Neveu-Schwarz (NS), and
periodic, also called Ramond (R), boundary conditions on the cycles of the
surface. The entropy of Jordan-Wigner fermions, instead, corresponds to
antiperiodic boundary conditions over all cycles (see also Ref. Camargo2021 ).
As discussed in Ref. Headrick2010  in detail, indeed, the discrepancy between
the entanglement entropy of two disjoint intervals for the free Dirac fermion
Casini2005 ; Casini2009reduced_density ; Casini2009Entanglement  and for the
free boson compactified on a circle with the corresponding radius
Calabrese2009Entanglement1 , in our notations $R=1$, can be traced back to the
spin structure. The result for the boson matches, instead, that for the
modular invariant Dirac fermion Calabrese2011Entanglement ; Coser2016Spin ,
which is obtained by summing over all boundary conditions. We mention that
similar problems related to Bose-Fermi duality and modular invariance for the
Dirac fermion on the torus arise when computing the entanglement of a single
interval at finite temperature Lokhande2015 ; Mukhi2018 ; Arias2020 .

The tripartite information in the XXZ model can be deduced from the results of
Refs Calabrese2009Entanglement1 ; Furukawa2009Mutual  on the free compactified
boson—see also Ref. Alba2011Entanglement . In light of Refs
Fagotti2010disjoint ; Coser2016Spin , we argue that it is possible to approach
this problem also in the fermionic picture by considering the massless
Thirring model. Concerning Ref. Coser2016Spin , we note that interacting spin
chains described by the free compactified boson with radius $R\neq 1$ have
been left out of the discussion. An analogous comment holds true for Ref.
Fagotti2010disjoint , which was focussed on noninteracting spin chains. Both
the construction of Ref. Fagotti2010disjoint  and the arguments exhibited in
Ref. Coser2016Spin  about the correspondence between partition functions and
moments of the pseudo-RDMs do not seem, however, to rely on the absence of
interactions ($R=1$). Thus we argue that the main observations of Refs
Fagotti2010disjoint ; Coser2016Spin  apply also to the Heisenberg XXZ chain,
whose low-energy properties are captured by the massless Thirring model, and,
in turn, by the free boson compactified on a circle. The problem of computing
the partition function of the massless Thirring model on the higher genus
Riemann surface with specified spin structure has been addressed in the works
by Freedman, Pilch Freedman1988Thirring ; Freedman1989Thirring , and Wu
Wu1989Determinants , by direct manipulations of fermionic path integrals. We
warn the reader that those results are not totally satisfactory: they display
inconsistencies about factors and signs, which were commented and partially
resolved in Ref. Sachs1996 . The tripartite information, however, can be
computed without knowing all the details of the partition functions and we
will fix the remaining ambiguities by imposing the known results for the XXZ
chain. These observations are the starting point of our investigations.

The first part of this paper is mainly derivative of Refs Freedman1989Thirring
; Calabrese2009Entanglement1 ; Fagotti2012New . On the one hand, we consider
the tripartite information in the XXZ fermion chain, which is described by the
massless Thirring model with antiperiodic boundary conditions on all cycles: A
direct application of Ref. Freedman1989Thirring  allows us to generalize the
analysis of Ref. Casini2005 , focussed on the free Dirac fermion, to the
interacting case, where we report a contrastingly nonzero tripartite
information. On the other hand, we reproduce the result for the XXZ spin
chain/free compactified boson obtained in Ref. Calabrese2009Entanglement1  by
summing over all the periodic and anti-periodic spin structures computed in
Ref. Freedman1989Thirring . We also consider a family of interacting spin-
chain models (with local translationally invariant Hamiltonians) with higher
central charge, that are related to XXZ by Jordan-Wigner transformations and
restrictions to sub-chains. We call the latter models Jordan-Wigner
deformations of XXZ. They are interacting generalizations of some free models
introduced in Ref. Suzuki1971The , whose tripartite information was studied in
Ref. Fagotti2012New . Their low energy description is a sum of independent
Thirring models with different species of fermions and correlated spin
structures. Knowing the contribution of each spin structure, we can easily
work out their tripartite information.

In the second part of the paper we set up the inverse problem of computing the
contribution of each spin structure from the knowledge of the tripartite
information in a family of spin chains. The simplest and most accessible case
is the torus, for which we can express the partition functions of the massless
Thirring model with given boundary conditions in terms of the Rényi-$2$
tripartite information of four JW deformations of XXZ. This yields a procedure
to check the QFT predictions against numerical simulations (e.g., with tensor
network algorithms). While effective on paper, this method is affected by a
problem that makes the numerical investigations challenging: both the range of
the interactions and the central charge of the models to consider increase by
$1$ for each additional deformation. That is to say, even restricting to the
torus, the numerical algorithm should be good enough to allow for the analysis
of the entropies of large spin blocks also in a model with range $4$ and
central charge $4$. The latter requirement is particularly challenging because
the entropy is proportional to the central charge, so the numerical
simulations become less and less effective Pollmann2009Theory ;
Stojevic2015Conformal . In order to circumvent those issues, we propose an
alternative procedure based on considering Hamiltonians obtained by
“interleaving” XXZ critical Hamiltonians with noncritical ones. This breaks
shift invariance from one-site into two-sites but provides a bound both on the
range of the interactions, which can be as small as $2$ for every model to
consider, and on their central charge, which remains the same as in the
original system (in the specific case of XXZ, it remains equal to $1$).
Remarkably, even if we focus on the torus and on the XXZ model, the procedure
that we propose can be applied to Riemann surfaces with higher genus (though,
it could become trickier) and to other spin-chain models with spin-flip
symmetry.

In the third part of the paper we go out of equilibrium and speculate about
the behaviour of the tripartite information at late times after having
prepared the system in a state that is close to but still different from the
ground state. Using the terminology of Ref. Calabrese2012Quantum , that would
correspond to a “small global quench” from another critical Hamiltonian. We
have recently investigated the tripartite information after global quenches in
quantum spin chains that can be mapped to free fermions Maric2022Universality
; Maric2023Universality ; Maric2023universality2 . In critical systems, in
particular, we showed that it can remain nonzero also in the limit of infinite
time and infinite lengths (taking the latter limit after the former one).
Remarkably, the tripartite information remains a function of the cross ratio
(1.4), but the symmetry $x\leftrightarrow 1-x$, which is always present in the
ground state, doesn’t survive the limits. Refs Maric2022Universality ;
Maric2023Universality ; Maric2023universality2  also showed that the limit
$x\rightarrow 1^{-}$ (taken after the limits of infinite time and lengths),
which was dubbed “residual tripartite information”, can be nonzero. This can
be surprising because it corresponds to the configuration in which the
distance between the blocks is negligible with respect to the blocks’ lengths,
therefore it could be naively identified with the degenerate case in which the
blocks are adjacent, for which the tripartite information is identically zero.
Clearly the limits do not commute, but there is something more: in all the
(noninteracting) cases investigated in Refs Maric2022Universality ;
Maric2023Universality ; Maric2023universality2 , the residual tripartite
information was either $0$ or $-\log 2$. Here we provide evidence that this is
not just a special feature of noninteracting models, we find indeed that also
the residual tripartite information after quenches of the anistropy in the
gapless phase of XXZ should be expected to be $-\log 2$.

##  2 Summary of results

We consider the XXZ spin-$\frac{1}{2}$ chain described by the Hamiltonian

$\bm{\mathbf{H}}_{\mathrm{XXZ}}=\sum_{\ell}\left(\bm{\mathbf{\sigma}}_{\ell}^{x}\bm{\mathbf{\sigma}}_{\ell+1}^{x}+\bm{\mathbf{\sigma}}_{\ell}^{y}\bm{\mathbf{\sigma}}_{\ell+1}^{y}+\Delta\bm{\mathbf{\sigma}}_{\ell}^{z}\bm{\mathbf{\sigma}}_{\ell+1}^{z}-2h\bm{\mathbf{\sigma}}_{\ell}^{z}\right)\,,$
(2.1)

where $h$ represents the coupling with an external magnetic field in the $z$
direction and $\Delta$ is usually referred to as “anisotropy”. In a strip of
the parameter space including $h=0\wedge|\Delta|<1$, the low-energy properties
of the model are described by a conformal field theory with central charge
$c=1$. By exploiting the correspondence between the XXZ spin-$\frac{1}{2}$
chain and the massless Thirring model, we point out the spin structures of the
former. This observation allows us to compute the tripartite information of
several related translationally invariant models. We report, in particular,
the results for the following models (see Figure 1 for a summary of a relation
between the tripartite information and Thirring partition functions in the
case $\alpha=2$)

  * •

XXZ fermion chain. This is a model of spinless interacting fermions, dual to
the XXZ spin chain by the Jordan-Wigner transformation
$\bm{\mathbf{c}}^{\dagger}_{\ell}=\prod_{j<\ell}\bm{\mathbf{\sigma}}_{j}^{z}\,\frac{\bm{\mathbf{\sigma}}_{\ell}^{x}+i\bm{\mathbf{\sigma}}_{\ell}^{y}}{2}$,
where $\bm{\mathbf{c}}^{\dagger}$ are creation operators of spinless fermions
satisfying
$\\{\bm{\mathbf{c}}^{\dagger}_{\ell},\bm{\mathbf{c}}_{n}\\}=\delta_{\ell
n}\bm{\mathbf{1}}$ and $\\{\bm{\mathbf{c}}_{\ell},\bm{\mathbf{c}}_{n}\\}=0$.
Its Hamiltonian reads

$\bm{\mathbf{H}}_{\mathrm{fXXZ}}=2\sum_{\ell}\left(\bm{\mathbf{c}}_{\ell+1}\bm{\mathbf{c}}^{\dagger}_{\ell}+\bm{\mathbf{c}}_{\ell}\bm{\mathbf{c}}^{\dagger}_{\ell+1}-2(\Delta+h)\bm{\mathbf{c}}^{\dagger}_{\ell}\bm{\mathbf{c}}_{\ell}+2\Delta\bm{\mathbf{c}}^{\dagger}_{\ell}\bm{\mathbf{c}}_{\ell}\bm{\mathbf{c}}^{\dagger}_{\ell+1}\bm{\mathbf{c}}_{\ell+1}\right)\,.$
(2.2)

  * •

Jordan-Wigner $\nu$-deformation of the XXZ spin chain. This is dual to $\nu$
uncoupled XXZ fermion sub-chains interleaved with one another

$\begin{split}\bm{\mathbf{H}}_{\mathrm{XXZ}}^{(\nu)}=&\sum_{\ell}\Big{[}\bm{\mathbf{\sigma}}_{\ell}^{x}(\bm{\mathbf{\sigma}}_{\ell+1}^{z}\bm{\mathbf{\sigma}}_{\ell+2}^{z}\cdots\bm{\mathbf{\sigma}}_{\ell+\nu-1}^{z})\bm{\mathbf{\sigma}}_{\ell+\nu}^{x}+\bm{\mathbf{\sigma}}_{\ell}^{y}(\bm{\mathbf{\sigma}}_{\ell+1}^{z}\bm{\mathbf{\sigma}}_{\ell+2}^{z}\cdots\bm{\mathbf{\sigma}}_{\ell+\nu-1}^{z})\bm{\mathbf{\sigma}}_{\ell+\nu}^{y}\\\
&+\Delta\bm{\mathbf{\sigma}}_{\ell}^{z}\bm{\mathbf{\sigma}}_{\ell+\nu}^{z}-2h\bm{\mathbf{\sigma}}_{\ell}^{z}\Big{]}\,.\end{split}$
(2.3)

$a$$b$ Torus with modular parameter $\tau=i\frac{K(1-x)}{K(x)}$
$Z_{\varepsilon,\delta}^{\mathrm{QFT}}\equiv$ partition function of fermionic
QFT on the torus, with boundary conditions $\varepsilon,\delta$ on cycle $a,b$
respectively Jordan-Wigner transformation Fermionic XXZGround
state$I_{3}^{(\alpha=2)}=\log\frac{Z_{\mathrm{NS},\mathrm{NS}}^{\mathrm{Thirring}}}{Z_{\mathrm{NS},\mathrm{NS}}^{\mathrm{Dirac}}}$
Small quantum quench Stationary
state$I_{3}^{(\alpha=2)}=\log\frac{Z_{\mathrm{NS},\mathrm{NS}}^{\mathrm{Thirring}}}{Z_{\mathrm{NS},\mathrm{NS}}^{\mathrm{Dirac}}}$XXZGround
state$I_{3}^{(\alpha=2)}=\log\sum\limits_{\varepsilon,\delta=\mathrm{NS},\mathrm{R}}\frac{(-1)^{4\epsilon\delta}Z_{\varepsilon,\delta}^{\mathrm{Thirring}}}{2Z_{\mathrm{NS},\mathrm{NS}}^{\mathrm{Dirac}}}$
Small quantum quench Stationary
state$I_{3}^{(\alpha=2)}=\log\sum\limits_{\delta=\mathrm{NS},\mathrm{R}}\frac{Z_{\mathrm{NS},\delta}^{\mathrm{Thirring}}}{2Z_{\mathrm{NS},\mathrm{NS}}^{\mathrm{Dirac}}}$
Figure 1: Summary of the relation between the tripartite information
$I_{3}^{(\alpha)}(x)$ and QFT partition functions for the Rényi index
$\alpha=2$. The Thirring coupling constant is fixed by the model parameters.
Figure 2: The Rényi-$2$ tripartite information $I_{3}^{(\alpha=2)}$ as a
function of the cross ratio $x$ in the ground state of the fermionic XXZ chain
for $h=0$ and different values of $\Delta$, reported in the legend. Depending
on the model parameters, the tripartite information can be positive, negative
or zero (the free case is $\Delta=0$).

XXZ fermion chain:

Contrary to free-fermion systems, such as the Kitaev model, the Rényi-$\alpha$
tripartite information in the XXZ fermion chain is not identically zero.
Specifically, we find

$I_{3}^{(\alpha)}(x)=\frac{1}{\alpha-1}\log\left(\sum_{\vec{\varepsilon}\in\\{0,\frac{1}{2}\\}^{\alpha-1}}\frac{\Theta\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{0}\end{smallmatrix}\right]\left(\vec{0}|\frac{1}{\eta}\Omega\right)\Theta\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{0}\end{smallmatrix}\right]\left(\vec{0}|4\eta\Omega\right)}{\left|\Theta\left(\vec{0}|\Omega\right)\right|^{2}}\right)\,.$
(2.4)

Here $\eta$ is the universal critical exponent proportional to the square of
the compactification radius of the underlying bosonic theory, which for $h=0$
reads $\eta=\frac{\arccos(-\Delta)}{\pi}$; $\Omega$ is the
$(\alpha-1)\times(\alpha-1)$ period matrix of the Riemann surface
$\mathcal{R}_{\alpha}$ with elements

$[\Omega]_{\ell n}=\frac{2i}{\alpha}\sum_{k=1}^{\alpha-1}\sin(\tfrac{\pi
k}{\alpha})\cos(\tfrac{2\pi
k(\ell-n)}{\alpha})\beta_{\frac{k}{\alpha}}(x),\quad\beta_{\frac{k}{\alpha}}(x)\equiv\tfrac{P_{-k/\alpha}(2x-1)}{P_{-k/\alpha}(1-2x)}\,,$
(2.5)

where $P_{\mu}(z)$ denotes the Legendre functions;
$\Theta[\vec{e}](\vec{z}|M)$ is the Riemann-Siegel theta function with
characteristic, i.e.,

$\Theta[\vec{e}](\vec{z}|M)=\sum_{\vec{m}\in\mathbb{Z}^{\alpha-1}}e^{i\pi(\vec{m}+\vec{\varepsilon})^{t}M(\vec{m}+\vec{\varepsilon})+2\pi
i(\vec{m}+\vec{\varepsilon})\cdot(\vec{z}+\vec{\delta})}\qquad\vec{e}=\begin{pmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{pmatrix}\,.$ (2.6)

Note that it is customary to omit the characteristic when it is zero, i.e.,
$\Theta\equiv\Theta[\vec{0}]$. The behavior of the tripartite information for
$\alpha=2$ is shown in Fig. 2. The curves for higher $\alpha$ are barely
distinguishable from those for $\alpha=2$ so we do not present them.

We also obtain the von Neumann tripartite information in the limit of small
$x$

$\begin{split}I_{3}^{(\mathrm{vN})}(x)=&-s_{\mathrm{vN}}(1)\frac{x}{2}+s_{\mathrm{vN}}(\tfrac{1}{\eta})\left(\frac{x}{4}\right)^{\tfrac{1}{\eta}}+s_{\mathrm{vN}}(4\eta)\left(\frac{x}{4}\right)^{4\eta}\\\
&+s_{\mathrm{vN}}(\tfrac{1}{4\eta}+\eta)\left(\frac{x}{4}\right)^{\tfrac{1}{4\eta}+\eta}+o\big{(}x^{\min(1,\tfrac{1}{\eta},4\eta)}\big{)},\qquad\textrm{as
}x\to 0,\end{split}$ (2.7)

where $s_{\mathrm{vN}}$ is defined as Calabrese2011Entanglement

$s_{\mathrm{vN}}(\mathcal{A})=\frac{\sqrt{\pi}\Gamma(\mathcal{A}+1)}{2\Gamma(\mathcal{A}+\frac{3}{2})},\quad\mathcal{A}>0.$
(2.8)

We find, in particular, that the tripartite information is negative for
$\frac{1}{4}\leq\eta\leq 1$ (i.e., $\Delta\geq-\frac{\sqrt{2}}{2}$ when
$h=0$), while in the XXZ spin chain it is always non-negative. For the
Rényi-$\alpha$ tripartite information with $\alpha=2,3$ we report an unusual
behavior: for certain values of $\Delta$, tripartite information changes the
sign as a function of $x$ (see Fig. 3).

Figure 3: The Rényi-$\alpha$ tripartite information $I_{3}^{(\alpha)}$ as a
function of the cross ratio $x$ in the ground state of the fermionic XXZ chain
for $\Delta=-0.76$ and $h=0$. The values of the Rényi index $\alpha$ are
reported in the legend. For $\alpha=2,3$ tripartite information changes sign
as a function of $x$.

(a) Contribution from different spin structures for $\alpha=2$ and
$(\Delta,h)=(-0.3,0)$. The three subsystems have equal lengths $|A|=|B|=|C|$,
so the cross ratio is $x=1/4$. The dashed lines show the free case
($\Delta=0$) for comparison.

(b) Contribution from different spin structures for $\alpha=2$ and
$(\Delta,h)=(-0.3,0)$. The subsystem lengths are set to be
$(|A|,|B|,|C|)=(2\ell,\ell,3\ell)$ for $\ell=1,2,3,\ldots$, so the cross ratio
is $x=1/2$. The lines R-NS and NS-R overlap. The dashed lines show the free
case for comparison.

Figure 4: Contributions of different spin structures $(\varepsilon,\delta)$ to
the Rényi-$2$ tripartite information of three finite adjacent subsystems
$A,B,C$ in the XXZ spin chain on infinite lattice. Namely, function
$\mathcal{F}_{f}\left[\begin{smallmatrix}\varepsilon\\\
\delta\end{smallmatrix}\right](A,B,C)$ (see Eq. (3.18) and Eq. (3.19)) is
shown for different subsystem lengths and all four torus spin structures. In
each plot the lengths of the three subsystems are chosen so that the cross
ratio $x=|A||C|/[(|A|+|B|)(|B|+|C|)]$ is fixed. The numerical data
(represented by points) is obtained by considering the method of fusion of
models, discussed in the text, and by incorporating tensor networks
simulations. As the subsystems lengths get bigger, the data can be seen to
converge towards the prediction (solid lines), given by Eq. (2.10), obtained
from the massless Thirring QFT considerations.

Expansion in the spin structure.

The tripartite information in the ground state of a spin-flip invariant
spin-$\frac{1}{2}$ chain Hamiltonian can be expanded as follows

$I_{3}^{(\alpha)}(A,B,C)=\frac{1}{\alpha-1}\log\left(\frac{1}{2^{\alpha-1}}\sum_{\vec{\varepsilon},\vec{\delta}\in\\{0,\frac{1}{2}\\}^{\alpha-1}}(-1)^{4\vec{\varepsilon}\cdot\vec{\delta}}\mathcal{F}_{f}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](A,B,C)\right)$ (2.9)

where $\mathcal{F}_{f}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](A,B,C)$ describes the contribution from
the underlying fermionic system with given spin structure—see Sec. 3.1.

XXZ spin chain.

In the Heisenberg XXZ model we find

$\mathcal{F}_{f}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](A,B,C)=\sum_{\vec{\mu}\in\\{0,\frac{1}{2}\\}^{\alpha-1}}(-1)^{4\vec{\mu}\cdot\vec{\delta}}\frac{\Theta\left[\begin{smallmatrix}\vec{\epsilon}+\vec{\mu}\\\
\vec{0}\end{smallmatrix}\right]\left(\vec{0}|\frac{1}{\eta}\Omega\right)\Theta\left[\begin{smallmatrix}\vec{\mu}\\\
\vec{0}\end{smallmatrix}\right]\left(\vec{0}|4\eta\Omega\right)}{\left[\Theta(\vec{0}|\Omega)\right]^{2}}\,.$
(2.10)

Note, in particular, that the fermionic tripartite information—(2.4)—is
obtained by taking only the term with $\vec{\varepsilon}=\vec{\delta}=0$.
While the sum over all spin structures has the well-known invariance under
$\eta\rightarrow\frac{1}{\eta}$ (this is shown explicitly in Appendix B), each
individual term exhibits another symmetry:

$\mathcal{F}_{f}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](A,B,C)\xrightarrow{\eta\rightarrow\frac{1}{4\eta}}(-1)^{4\vec{\varepsilon}\cdot\vec{\delta}}\mathcal{F}_{f}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](A,B,C)\,.$ (2.11)

We point out that the transformation $\eta\rightarrow\frac{1}{4\eta}$
distinguishes the even spin structures from the odd ones. Here we call a spin
structure even or odd depending on whether
$(-1)^{4\vec{\varepsilon}\cdot\vec{\delta}}=1$ or
$(-1)^{4\vec{\varepsilon}\cdot\vec{\delta}}=-1$ respectively. For $h=0$ this
transformation connects two XXZ chains as long as
$\Delta\geq-\frac{\sqrt{2}}{2}$. For $h=0$ and $\Delta<-\frac{\sqrt{2}}{2}$ a
similar effect can be achieved by applying a double transformation
$\eta\rightarrow\frac{1}{4\eta}\rightarrow 4\eta$, which doesn’t preserve the
terms separately but, effectively, acts as it reversed the sign of the odd
spin structures.

Jordan-Wigner deformations:

We work out the tripartite information in the Jordan-Wigner
deformations—(2.3)—of the XXZ spin chain. In agreement with the discussion of
Ref. Fagotti2012New , we find that the tripartite information is given by

$I_{3}^{(\alpha;\nu)}(x)=\frac{1}{\alpha-1}\log\left(\frac{1}{2^{\alpha-1}}\sum_{\vec{\varepsilon},\vec{\delta}\in\\{0,\frac{1}{2}\\}^{\alpha-1}}(-1)^{4\vec{\varepsilon}\cdot\vec{\delta}}\mathcal{F}_{f}^{(\nu)}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](x)\right)$ (2.12)

with

$\mathcal{F}_{f}^{(\nu)}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](x)=\Bigl{(}\mathcal{F}_{f}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](x)\Bigr{)}^{\nu}\,,$ (2.13)

where $\mathcal{F}_{f}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](x)$ is given by Eq. (2.10). For $\nu=2$
the spin tripartite information matches the fermionic one, which, in turn, is
equal to twice the tripartite information in the fermionic XXZ chain

$I_{3}^{(\alpha;\nu=2)}(x)=2\ I_{3}^{(\alpha;\mathrm{fXXZ})}(x)\;,$ (2.14)

where $I_{3}^{(\alpha;\mathrm{fXXZ})}$ is given by Eq. (2.4). For small $x$
and any $\nu$ we have computed the von Neumann tripatite information. The
results are given in Sec. 5.2. The tripartite information can be positive,
negative or zero, depending on the model parameters.

Figure 5: Prediction for the Rényi-$2$ tripartite information
$I_{3}^{(\alpha=2)}$ as a function of the cross ratio $x$ in the stationary
state after the small quantum quench in the XXZ spin chain for $h=0$ and
different values of $\Delta$, reported in the legend. The tripartite
information depends on the model parameters, but, irrespectively of them, it
reaches the value $-\log 2$ as $x\to 1^{-}$.

From tripartite information to spin structures:

We propose two different procedures for computing the contribution of each
spin structure of the massless Thirring model on a Riemann surface to the
Rényi-$\alpha$ tripartite information in the XXZ chain, with the main focus on
the torus, i.e. $\alpha=2$.

Jordan-Wigner $\nu$-deformation of the XXZ spin chain We express the partition
function with given spin structure in terms of the Rényi-$2$ tripartite
information of four JW deformations of XXZ (including the XXZ model itself).
For example, the term corresponding to Ramond boundary conditions on both
cycles, which we have indicated by
$\mathcal{F}_{f}\left[\begin{smallmatrix}\nicefrac{{1}}{{2}}\\\
\nicefrac{{1}}{{2}}\end{smallmatrix}\right](x)$, can be written as follows

$\mathcal{F}_{f}\left[\begin{smallmatrix}\nicefrac{{1}}{{2}}\\\
\nicefrac{{1}}{{2}}\end{smallmatrix}\right](x)=\tfrac{4e^{4I_{3}^{(2;1)}(x)}-12e^{2I_{3}^{(2;1)}(x)}e^{I_{3}^{(2;2)}(x)}+3e^{2I_{3}^{(2;2)}(x)}+8e^{I_{3}^{(2;1)}(x)}e^{I_{3}^{(2;3)}(x)}-3e^{I_{3}^{(2;4)}(x)}}{12e^{I_{3}^{(2;1)}(x)}e^{I_{3}^{(2;2)}(x)}-8e^{3I_{3}^{(2;1)}(x)}-4e^{I_{3}^{(2;3)}(x)}}\,,$
(2.15)

where $I_{3}^{(2;\nu)}(x)$ denotes the Rényi-2 tripartite information of the
JW $\nu$-deformation of XXZ. In fact, a posteriori we observe that this
particular term can be obtained just studying the tripartite information in
XXZ chains. By exploiting the symmetry mentioned before, we indeed find

$\mathcal{F}_{f}^{(\eta)}\left[\begin{smallmatrix}\nicefrac{{1}}{{2}}\\\
\nicefrac{{1}}{{2}}\end{smallmatrix}\right](x)=e^{I_{3}^{(2)}(x;\eta)}-e^{I_{3}^{(2)}(x;\min(\nicefrac{{1}}{{(4\eta)}},4\eta))}\,,$
(2.16)

where the subscripts $\eta$, $4\eta$, and $\nicefrac{{1}}{{(4\eta)}}$ identify
the XXZ model.

Interleaved models The second procedure is based on the study of a class of
two-site shift invariant Hamiltonians such as the following

$\bm{\mathbf{H}}=\bm{\mathbf{H}}_{e}+\bm{\mathbf{H}}_{o}$ (2.17)

where $\bm{\mathbf{H}}_{e}$ and $\bm{\mathbf{H}}_{o}$ are commuting
Hamiltonians, $[\bm{\mathbf{H}}_{e},\bm{\mathbf{H}}_{o}]=0$, with

$\displaystyle\bm{\mathbf{H}}_{e}=$
$\displaystyle\sum_{\ell}\left(\bm{\mathbf{\sigma}}_{2\ell}^{x}\bm{\mathbf{\sigma}}_{2\ell+1}^{z}\bm{\mathbf{\sigma}}_{2\ell+1}^{x}+\bm{\mathbf{\sigma}}_{2\ell}^{y}\bm{\mathbf{\sigma}}_{2\ell+1}^{z}\bm{\mathbf{\sigma}}_{2\ell+1}^{y}+\Delta\bm{\mathbf{\sigma}}_{2\ell}^{z}\bm{\mathbf{\sigma}}_{2\ell+2}^{z}-2h\bm{\mathbf{\sigma}}_{2\ell}^{z}\right)\,$
(2.18) $\displaystyle\bm{\mathbf{H}}_{o}=$
$\displaystyle-\sum_{\ell}\left[(1+\gamma)\bm{\mathbf{\sigma}}_{2\ell-1}^{x}\bm{\mathbf{\sigma}}_{2\ell}^{z}\bm{\mathbf{\sigma}}_{2\ell+1}^{x}+(1-\gamma)\bm{\mathbf{\sigma}}_{2\ell-1}^{y}\bm{\mathbf{\sigma}}_{2\ell}^{z}\bm{\mathbf{\sigma}}_{2\ell+1}^{y}-2h^{\prime}\bm{\mathbf{\sigma}}_{2\ell-1}^{z}\right]\,.$
(2.19)

This describes an XXZ model interleaved with an XY model, and the ground state
is the tensor product of two Fock states associated with independent sets of
Jordan-Wigner fermions.

We point out

$\mathcal{F}_{f}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](A,B,C)=\mathcal{F}_{f}^{(\mathrm{XXZ}_{\Delta,h})}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](A_{e},B_{e},C_{e})\mathcal{F}_{f}^{(\mathrm{XY}_{\gamma,h^{\prime}})}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](A_{o},B_{o},C_{o})\,,$ (2.20)

where we used the notation $X_{e(o)}$ to indicate the set of even (odd) sites
of subsystem $X$ and the elements on the right hand side correspond to the
underlying XXZ—(2.10)—and XY chains— Section 3.3. As a result, the ground
state tripartite information can still be written as a linear combination of
Thirring partition functions with given spin structure, the coefficients of
the combination depending on $\gamma$ and $h^{\prime}$. We exploit this
connection with the Thirring partition functions to perform a numerical check
of the analytical predictions for given spin structure. Some results are
presented in Fig. 4. To the best of our knowledge, this is the first numerical
analysis of its kind in an interacting spin chain (we note, in particular,
that the numerical analysis of Refs Fagotti2012New ; Coser2016Spin  had the
same goal but the technique that was used is specific to noninteracting spin
chains).

Global quenches from critical states:

Refs Maric2022Universality ; Maric2023Universality  studied the stationary
value of the tripartite information long after having prepared the system in a
critical state (a collection of Fermi seas). It was in particular observed
that the main effect of the quench is to turn some algebraically decaying
correlations into exponentially decaying ones. The correlations that become
exponentially small are those associated with the vertex operator in the
underlying free boson compactified on a circle. Assuming this to remain true
also in the interacting case, we speculate about the behaviour of the
tripartite information after global quenches. This corresponds to dropping all
spin structures with $\vec{\varepsilon}\neq\vec{0}$ and we conjecture, in
particular, a formula for the tripartite information in the generalised Gibbs
ensemble emerging after a quench of the anisotropy
$\Delta_{0}\rightarrow\Delta$ in the XXZ spin-$\frac{1}{2}$ chain in the limit
of small quench

$\lim_{\Delta_{0}\rightarrow\Delta}\lim_{t\rightarrow\infty}I_{3}^{(\alpha)}(x)=\frac{1}{\alpha-1}\log\frac{\Theta(\vec{0}|\frac{1}{\eta}\Omega)\Theta(\vec{0}|4\eta\Omega)}{[\Theta(\vec{0}|\Omega)]^{2}}\,.$
(2.21)

The behavior for $\alpha=2$ is presented in Fig. 5. We have not been able to
compare this conjecture with data from numerical simulations, but the interest
aroused by nonequilibrium time evolution in lattice systems is triggering
rapid development of efficient algorithms that could make a numerical
verification possible before long.

From Eq. (2.21) we have been able to obtain the large $x$ expansion of the von
Neumann tripartite information,

$\begin{split}&\lim_{\Delta_{0}\rightarrow\Delta}\lim_{t\rightarrow\infty}I_{3}^{(\mathrm{vN})}(x)\\\
&\xrightarrow{1-x\ll 1}-\log 2\
+s_{\mathrm{vN}}(\eta)\left(\frac{1-x}{4}\right)^{\eta}+s_{\mathrm{vN}}(\tfrac{1}{4\eta})\left(\frac{1-x}{4}\right)^{\frac{1}{4\eta}}+o\left((1-x)^{\min(\eta,\tfrac{1}{4\eta})}\right),\end{split}$
(2.22)

where the function $s_{\mathrm{vN}}$ is defined in Eq. (2.8).

Eqs (2.21) and (2.22) allow us to compute the residual tripartite information,
which is the limit $x\rightarrow 1^{-}$. Remarkably, we find that it is
independent of the anisotropy and equals $-\log 2$. Even if our analysis is
based on a conjecture that deserves further investigation, our findings
support the claim of universality made in Ref. Maric2022Universality , which
was based on the study of a wide class of noninteracting spin-chain models.

##  3 Tripartite information

###  3.1 Fermions vs spins: spin structures

In spin-lattice systems, the tensor product structure of the Hilbert space
allows for a clear notion of separable states and entanglement. In composite
fermionic systems the Fock space does not exhibit the same structure, and the
situation is not as simple. We refer the reader to Ref. Banuls2007  for a
discussion on this topic. In this section we explore the differences between
the entanglement entropies of blocks of spins in spin-$\frac{1}{2}$ chains
with spin-flip symmetry and the entropies of the corresponding spinless
fermions defined through a Jordan-Wigner transformation.

#### Reduced density matrix.

In a spin-$\frac{1}{2}$ system, such as the XXZ chain, the sites of a
subsystem $X$ are associated with the indices of the spin operators
$\bm{\mathbf{\sigma}}_{\ell}^{\gamma}$, $\ell\in X$, for $\gamma=x,y,z$. If
the whole system is in a state $\ket{\Psi}$, the reduced density matrix
$\rho_{X}$ of a subsystem $X$ reads

$\rho_{X}=\frac{1}{2^{|X|}}\sum_{\gamma_{\ell}\in\\{0,x,y,z\\},\ell\in
X}\bra{\Psi}\prod_{\ell\in
X}\bm{\mathbf{\sigma}}^{\gamma_{\ell}}_{\ell}\ket{\Psi}\bigotimes_{\ell\in
X}\sigma^{\gamma_{\ell}}\;,$ (3.1)

where the sum is over an orthogonal basis of operators on $X$, given by all
possible products of Pauli matrices (referred to in the following as
“strings”), and we use the standard convention $\sigma^{0}\equiv\mathbb{I}$.
Note that we use bold fonts for operators defined in the Hilbert space of the
whole lattice and unbold fonts for the ones defined in the Hilbert space
associated to a subsystem or a single site.

In systems of spinless fermions one can associate a subsystem
$X=\\{\ell_{1},\ell_{2},\ldots,\ell_{|X|}\\}$ with the indices of the
fermionic operators $\bm{\mathbf{c}}_{\ell}^{\dagger},\bm{\mathbf{c}}_{\ell}$
that create and annihilate, respectively, a fermion at site $\ell$. In
general, such a notion of a fermionic subsystem is problematic: the operators
constructed with fermions in the set $\\{\bm{\mathbf{c}}_{\ell}|\ell\in X\\}$
do not commute with every operator in the complementary set
$\\{\bm{\mathbf{c}}_{\ell}|\ell\in\bar{X}\\}$. However, we deal with states in
which expectation values of products of an odd number of fermionic operators
vanish and such problems do not arise. By introducing the fermionic operators
$c_{\ell}^{\dagger},c_{\ell}$, $\ell\in X$, that are defined on the $2^{|X|}$
dimensional Fock space asscociated to the subsystem (note the unbold
notation), one can define the operator

$\rho_{f,X}=\frac{1}{2^{|X|}}\sum_{F_{X}}\braket{\Psi}{\bm{\mathbf{F}}_{X}}{\Psi}F_{X}^{\dagger}\,,$
(3.2)

where we have denoted

$\bm{\mathbf{F}}_{X}=(\bm{\mathbf{c}}_{\ell_{1}}^{\dagger})^{\gamma_{1}}(\bm{\mathbf{c}}_{\ell_{1}})^{\delta_{1}}\ldots(\bm{\mathbf{c}}_{\ell_{|X|}}^{\dagger})^{\gamma_{{|X|}}}(\bm{\mathbf{c}}_{\ell_{|X|}})^{\delta_{{|X|}}}\;,\quad
F_{X}=(c_{\ell_{1}}^{\dagger})^{\gamma_{1}}(c_{\ell_{1}})^{\delta_{1}}\ldots(c_{\ell_{|X|}}^{\dagger})^{\gamma_{{|X|}}}(c_{\ell_{|X|}})^{\delta_{{|X|}}}\;,$
(3.3)

with $\gamma_{\ell},\delta_{\ell}\in\\{0,1\\},\ell\in X$ a generic fermionic
operator on the subsystem. The sum in (3.2) is over all such operators, i.e.
over all $\gamma_{\ell},\delta_{\ell}\in\\{0,1\\},\ell\in X$. Since we deal
with states where the expectation value of any product of an odd number of
fermions vanishes, the sum can be further restricted with the condition
$(-1)^{\sum_{j}\gamma_{j}}=(-1)^{\sum_{j}\delta_{j}}$. The operator (3.2) is a
fermionic reduced density matrix in the sense that it reproduces the
expectation values of operators associated with the subsystem,

$\mathrm{tr}\left(\rho_{f,X}F_{X}\right)=\bra{\Psi}\bm{\mathbf{F}}_{X}\ket{\Psi}.$
(3.4)

The entanglement entropies of a subsystem $X$ of fermions can then be computed
in an analogous way to spins, through
$S_{\alpha}(X)=\log[\mathrm{tr}(\rho_{f,X}^{\alpha})]/(1-\alpha)$ and
$S_{1}(X)=-\mathrm{tr}[\rho_{f,X}\log\rho_{f,X}]$. One can also take it as a
definition of fermionic entanglement entropies.

#### Jordan-Wigner transformation.

In view of comparing the entropies of spins and fermions under the Jordan-
Wigner transformation, due to its non-local properties it is more convenient
to work with the Hilbert space associated to the whole lattice. Thus, instead
of the reduced density matrices, we introduce the operators,

$\begin{split}\bm{\mathbf{R}}_{X}(\mathrm{spins})&=\mathrm{tr}[\bm{\mathbf{1}}_{X}]\rho_{X}\otimes
I_{\bar{X}}\\\ &=\sum_{\gamma_{\ell}\in\\{0,x,y,z\\},\ell\in
X}\braket{\Psi}{\prod_{\ell\in
X}\bm{\mathbf{\sigma}}^{\gamma_{\ell}}_{\ell}}{\Psi}\prod_{\ell\in
X}\bm{\mathbf{\sigma}}_{\ell}^{\gamma_{\ell}}\end{split}$ (3.5)

for spins, and

$\begin{split}\bm{\mathbf{R}}_{X}&=\sum_{\bm{\mathbf{F}}_{X}}\braket{\Psi}{\bm{\mathbf{F}}_{X}}{\Psi}\bm{\mathbf{F}}_{X}^{\dagger}\\\
&=\\!\\!\\!\\!\\!\\!\sum_{\gamma_{\ell},\delta_{\ell}\in\\{0,1\\},\ell\in
X\atop(-1)^{\sum_{j}\gamma_{j}}=(-1)^{\sum_{j}\delta_{j}}}\\!\\!\\!\\!\\!\\!\braket{\Psi}{(\bm{\mathbf{c}}_{\ell_{1}}^{\dagger})^{\gamma_{1}}(\bm{\mathbf{c}}_{\ell_{1}})^{\delta_{1}}\ldots(\bm{\mathbf{c}}_{\ell{X|}}^{\dagger})^{\gamma_{{|X|}}}(\bm{\mathbf{c}}_{\ell_{|X|}})^{\delta_{{|X|}}}|\Psi}{(}\bm{\mathbf{c}}^{\dagger}_{\ell_{|X|}})^{\delta_{{|X|}}}(\bm{\mathbf{c}}_{\ell_{|X|}})^{\gamma_{{|X|}}}\cdots(\bm{\mathbf{c}}_{\ell_{1}}^{\dagger})^{\delta_{1}}(\bm{\mathbf{c}}_{\ell_{1}})^{\gamma_{1}}\end{split}$
(3.6)

for fermions. Note that (3.1) and (3.5) do not differ only in the
normalisation, but also in the space where they are defined, as we emphasized
by using or not using bold fonts, and similarly for (3.2) and (3.6). The
entanglement entropies in terms of $\bm{\mathbf{R}}_{X}$ read

$\displaystyle S_{\alpha}(X)=$
$\displaystyle\frac{1}{1-\alpha}\log\frac{\mathrm{tr}[\bm{\mathbf{R}}_{X}^{\alpha}]}{\mathrm{tr}[\bm{\mathbf{1}}]}+\log\mathrm{tr}[\bm{\mathbf{1}}_{X}]\equiv\frac{1}{1-\alpha}\log\braket{\Psi}{\bm{\mathbf{R}}_{X}^{\alpha-1}}{\Psi}+\log\mathrm{tr}[\bm{\mathbf{1}}_{X}]$
(3.7) $\displaystyle S_{1}(X)=$
$\displaystyle-\frac{\mathrm{tr}[\bm{\mathbf{R}}_{X}\log\bm{\mathbf{R}}_{X}]}{\mathrm{tr}[\bm{\mathbf{1}}]}+\log\mathrm{tr}[\bm{\mathbf{1}}_{X}]\equiv-\braket{\Psi}{\log\bm{\mathbf{R}}_{X}}{\Psi}+\log\mathrm{tr}[\bm{\mathbf{1}}_{X}]\,.$

The Jordan-Wigner transformation
$\bm{\mathbf{\sigma}}_{\ell}^{+}=\prod_{j<\ell}(2\bm{\mathbf{c}}^{\dagger}_{j}\bm{\mathbf{c}}_{j}-1)\bm{\mathbf{c}}^{\dagger}_{\ell}$
maps the spins to spinless fermions and vice versa. The relation between the
entropy of the subsystem $X$ of spins and the corresponding one of the Jordan-
Wigner fermions depends on whether $X$ is connected or not.

  * •

If $X$ is a connected set of sites, the auxiliary operators for spins and
fermions are equal Vidal2003 ; Peschel2003 ; Peschel2009 ; Latorre2004 , i.e.
$\bm{\mathbf{R}}_{X}(\mathrm{spins})=\bm{\mathbf{R}}_{X}$. The Rényi entropies
of a single interval of spins and those of the corresponding interval of
Jordan-Wigner fermions are, in turn, equal. Note that if we use the Jordan-
Wigner transformation in the $2^{|X|}$ dimensional Hilbert space associated to
the subsystem $X$, we obtain

$\big{(}\bigotimes_{j=1}^{\ell-1}\sigma^{z}\big{)}\otimes\sigma^{+}\otimes\big{(}\bigotimes_{j=\ell+1}^{|X|}\mathbb{I}\big{)}=\prod_{j<\ell}(2c^{\dagger}_{j}c_{j}-1)c^{\dagger}_{\ell}\,,$

and hence we can adopt the point of view that the reduced density matrix of
spins is equal to the fermionic one, given by Eq. (3.2).

  * •

When $X$ is not connected, the non-locality of the Jordan-Wigner
transformation starts having effects. In particular, it has been observed in
Ref. Fagotti2010disjoint  that the Renyi-$\alpha$ entanglement entropy of two
disjoint intervals $X=A\cup C$ for a spin system in a state $\ket{\Psi}$ is
equal to the entropy of the corresponding Jordan-Wigner fermions in the state

$\ket{\Psi^{\prime}}=\left[\frac{1+(-1)^{\bm{\mathbf{F}}_{B}}}{2}+\frac{1-(-1)^{\bm{\mathbf{F}}_{B}}}{2}(-1)^{\bm{\mathbf{F}}_{C}}\right]\ket{\Psi},$
(3.8)

where $\bm{\mathbf{F}}_{B}=\sum_{\ell\in
B}\bm{\mathbf{c}}^{\dagger}_{\ell}\bm{\mathbf{c}}_{\ell}$ and
$\bm{\mathbf{F}}_{C}=\sum_{\ell\in
C}\bm{\mathbf{c}}^{\dagger}_{\ell}\bm{\mathbf{c}}_{\ell}$ is the fermion
number on interval $B$ and $C$ respectively and
$(-1)^{\bm{\mathbf{F}}}\equiv\exp(i\pi\bm{\mathbf{F}})$. That is to say, the
operator $\bm{\mathbf{R}}_{X}$ for spins, defined in (3.5), can be expressed
as a sum of four fermionic operators,

$\bm{\mathbf{R}}_{X}(\mathrm{spins})=\frac{\bm{\mathbf{R}}^{++}+\bm{\mathbf{R}}^{+-}}{2}+(-1)^{\bm{\mathbf{F}}_{B}}\frac{\bm{\mathbf{R}}^{-+}-\bm{\mathbf{R}}^{--}}{2},$
(3.9)

where

$\bm{\mathbf{R}}^{++}=\bm{\mathbf{R}}_{AC},\qquad\bm{\mathbf{R}}^{-+}=\sum_{F_{X}}\braket{\Psi}{(-1)^{\bm{\mathbf{F}}_{B}}\bm{\mathbf{F}}_{X}}{\Psi}F_{X}^{\dagger},\qquad\bm{\mathbf{R}}^{s-}=(-1)^{\bm{\mathbf{F}}_{C}}\bm{\mathbf{R}}^{s+}(-1)^{\bm{\mathbf{F}}_{C}},\,s=\pm
1,$ (3.10)

and $\bm{\mathbf{R}}_{AC}$ is defined in (3.6). Furthermore, the string
operator $(-1)^{\bm{\mathbf{F}}_{B}}$ in (3.9) can be removed without any
effect on the entanglement entropies. Namely, $\bm{\mathbf{R}}_{X}$ can be
replaced by the operator

$\bm{\mathbf{R}}^{\prime}_{X}(\mathrm{spins})=\frac{\bm{\mathbf{R}}^{++}+\bm{\mathbf{R}}^{+-}+\bm{\mathbf{R}}^{-+}-\bm{\mathbf{R}}^{--}}{2}$
(3.11)

in Eq. (3.7), indeed
$\mathrm{tr}[(\bm{\mathbf{R}}_{X}^{\prime}(\mathrm{spins}))^{\alpha}]=\mathrm{tr}[(\bm{\mathbf{R}}_{X}(\mathrm{spins})^{\alpha}]$.
This allows us to express the entanglement entropies in a less pedantic way
than through the auxiliary operators $\bm{\mathbf{R}}_{X}$, with the help of
the Jordan-Wigner transformation in the $2^{|A|+|C|}$ dimensional Hilbert
space associated with the subsystem $A\cup C$. The exact relations are
reported in the forthcoming Eq. (3.17).

#### Example.

We illustrate the difference between the entanglement of spins and fermions,
related through the Jordan-Wigner transformation, with a simple example. We
consider a system of four qubits, dual to a system of four fermions. We denote
the four adjacent sites by $A,B,C,D$ and consider the bipartition of the
system into subsystems $A\cup C$ and $B\cup D$, thus each consisting of two
disjoint sites.

The state

$\ket{\Psi_{1}}=\frac{1}{\sqrt{2}}\left(\ket{\downarrow\downarrow}+\ket{\uparrow\uparrow}\right)_{AC}\otimes\frac{1}{\sqrt{2}}\left(\ket{\downarrow\downarrow}+\ket{\uparrow\uparrow}\right)_{BD}$
(3.12)

is separable, i.e., it is a product state with respect to the bipartition into
$A\cup C$ and $B\cup D$ ($S_{\alpha}(A\cup C)=0$). In terms of fermions,
instead, the state reads

$\ket{\Psi_{1}}=\frac{1}{2}(1+\bm{\mathbf{c}}_{A}^{\dagger}\bm{\mathbf{c}}_{C}^{\dagger}+\bm{\mathbf{c}}_{B}^{\dagger}\bm{\mathbf{c}}_{D}^{\dagger}+\bm{\mathbf{c}}_{A}^{\dagger}\bm{\mathbf{c}}_{B}^{\dagger}\bm{\mathbf{c}}_{C}^{\dagger}\bm{\mathbf{c}}_{D}^{\dagger})\ket{0}.$
(3.13)

Roughly speaking, the state is not a fermionic product state, as witnessed by
the inequality
$\braket{(\bm{\mathbf{c}}_{A}\bm{\mathbf{c}}_{C})(\bm{\mathbf{c}}_{B}\bm{\mathbf{c}}_{D})}_{\Psi_{1}}\neq\braket{\bm{\mathbf{c}}_{A}\bm{\mathbf{c}}_{C}}_{\Psi_{1}}\braket{\bm{\mathbf{c}}_{B}\bm{\mathbf{c}}_{D}}_{\Psi_{1}}$,
in which only the left hand side is non-zero. The fermionic entropy
$S_{\alpha,f}(A\cup C)$ is indeed nonzero and equals $\log 2$.

On the other hand, the state

$\ket{\Psi_{2}}=\frac{1}{2}(1+\bm{\mathbf{c}}_{A}^{\dagger}\bm{\mathbf{c}}_{C}^{\dagger})(1+\bm{\mathbf{c}}_{B}^{\dagger}\bm{\mathbf{c}}_{D}^{\dagger})\ket{0}$
(3.14)

can be interpreted as a fermionic product state and has zero fermionic
entanglement entropy $S_{\alpha,f}(A\cup C)=0$. In terms of spins, it reads

$\ket{\Psi_{2}}=\frac{1}{\sqrt{2}}\left[\ket{\downarrow\downarrow}_{AC}\otimes\frac{1}{\sqrt{2}}(\ket{\downarrow\downarrow}+\ket{\uparrow\uparrow})_{BD}+\ket{\uparrow\uparrow}_{AC}\otimes\frac{1}{\sqrt{2}}(\ket{\downarrow\downarrow}-\ket{\uparrow\uparrow})_{BD}\right],$
(3.15)

which is instead an entangled state. The entropy of spins $S_{\alpha}(A\cup
C)$ now equals $\log 2$. One could easily see that relation (3.8) gives the
correct entropy in both examples.

#### Spin structures: lattice.

The spin RDM $\rho_{AC}$ (defined in (3.1)) is equivalent to the sum of four
pseudo-RDMs

$\rho_{AC}^{\prime}=\frac{\rho^{++}_{AC}+\rho_{AC}^{+-}+\rho_{AC}^{-+}-\rho_{AC}^{--}}{2}\,,$
(3.16)

in the sense that the two matrices have the same Rényi-$\alpha$ entanglement
entropies
($\mathrm{tr}[\rho_{AC}^{\alpha}]=\mathrm{tr}[(\rho_{AC}^{\prime})^{\alpha}]$).
The four pseudo-RDMs are defined as follows

$\begin{split}&\rho_{AC}^{+s^{\prime}}=\frac{1}{2^{|A|+|C|}}\left[\sum_{\mathrm{even}}\braket{\bm{\mathbf{O}}_{A}\bm{\mathbf{O}}_{C}}O_{A}\otimes
O_{C}+s^{\prime}\sum_{\mathrm{odd}}\braket{\bm{\mathbf{O}}_{A}(-1)^{\bm{\mathbf{F}}_{B}}\bm{\mathbf{O}}_{C}}O_{A}\otimes
O_{C}\right],\quad s^{\prime}=\pm 1\\\
&\rho_{AC}^{-s^{\prime}}=\frac{1}{2^{|A|+|C|}}\left[\sum_{\mathrm{even}}\braket{\bm{\mathbf{O}}_{A}(-1)^{\bm{\mathbf{F}}_{B}}\bm{\mathbf{O}}_{C}}O_{A}\otimes
O_{C}+s^{\prime}\sum_{\mathrm{odd}}\braket{\bm{\mathbf{O}}_{A}\bm{\mathbf{O}}_{C}}O_{A}\otimes
O_{C}\right],\quad s^{\prime}=\pm 1\end{split}$ (3.17)

where $O_{A}=\bigotimes_{\ell\in A}\sigma^{\gamma_{\ell}}$
($\bm{\mathbf{O}}_{A}=\prod_{\ell\in
A}\bm{\mathbf{\sigma}}_{\ell}^{\gamma_{\ell}}$ ) denotes a generic string of
Pauli operators on subsystem $A$, $O_{C}$ ($\bm{\mathbf{O}}_{C}$) denotes an
analogous string for subsystem $C$, $\sum_{\mathrm{even}}$ denotes the sum
over all $O_{A}$ and $O_{C}$ which are even under
$S^{z}_{A}\equiv\bigotimes_{\ell\in A}\sigma^{z}$ and
$S^{z}_{C}\equiv\bigotimes_{\ell\in C}\sigma^{z}$ respectively
($[O_{A},S_{A}^{z}]=0,[O_{C},S_{C}^{z}]=0$), and $\sum_{\mathrm{odd}}$ denotes
the sum over all $O_{A}$ and $O_{C}$ which are odd under $S^{z}_{A}$ and
$S^{z}_{C}$ respectively ($\\{O_{A},S_{A}^{z}\\}=0,\\{O_{C},S_{C}^{z}\\}=0$).

Note that both $\rho_{AC}^{++}$ and $\rho_{AC}^{+-}$ are proper RDMs, they are
indeed equivalent to the fermionic reduced density matrix (3.2). The matrices
$\rho_{AC}^{\pm s^{\prime}}$, $s^{\prime}=\pm 1$, on the other hand, cannot be
interpreted as density matrices because they are not positive semidefinite.
This is why, in general, we are referring to them as “pseudo-RDMs” rather than
RDMs.

The Rényi-$\alpha$ tripartite information can be written as follows:

$I_{3}^{(\alpha)}(A,B,C)=\frac{1}{\alpha-1}\log\left(\frac{1}{2^{\alpha-1}}\sum_{\vec{\varepsilon},\vec{\delta}\in\\{0,\frac{1}{2}\\}^{\alpha-1}}(-1)^{4\vec{\varepsilon}\cdot\vec{\delta}}\mathcal{F}_{f}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](A,B,C)\right)$ (3.18)

where

$\mathcal{F}_{f}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](A,B,C)=\frac{\mathrm{tr}[\rho_{AB}^{\alpha}]\mathrm{tr}[\rho_{BC}^{\alpha}]}{\mathrm{tr}[\rho_{A}^{\alpha}]\mathrm{tr}[\rho_{B}^{\alpha}]\mathrm{tr}[\rho_{C}^{\alpha}]\mathrm{tr}[\rho_{ABC}^{\alpha}]}\mathrm{tr}\left[\prod\limits_{k=1}^{\alpha}\rho_{AC}^{s_{k}[\vec{\varepsilon}],s_{k}^{\prime}[\vec{\delta}]}\right]$
(3.19)

with

$\displaystyle s_{k}[\vec{\varepsilon}]=$
$\displaystyle\begin{cases}(1-4\varepsilon_{1})\,,&k=1\\\
(1-4\varepsilon_{k})(1-4\varepsilon_{k-1})\,,&1<k<\alpha\\\
(1-4\varepsilon_{\alpha-1})\,,&k=\alpha\end{cases}$ (3.20) $\displaystyle
s^{\prime}_{k}[\vec{\delta}]=$ $\displaystyle\begin{cases}1\,,&k=1\\\
\prod\limits_{j=1}^{k-1}(1-4\delta_{j})\,,&1<k\leq\alpha\end{cases}\,.$

Note that, due to more or less explicit symmetries, there are alternative
equivalent identifications (e.g., we could change the sign of all
$s_{k}^{\prime}$). Incidentally, we isolated the sign
$(-1)^{4\vec{\varepsilon}\cdot\vec{\delta}}$ both because it naturally emerges
in the expansion in pseudo-RDMs and because it brings some advantages that
will be discussed later.

Eq. (3.18) is exact on the lattice. Ref. Fagotti2010disjoint  exploited this
structure to compute the Rényi entropies in spin systems that are mapped to
free fermions by the Jordan-Wigner transformation. Indeed, in that case,
$\rho^{s,s^{\prime}}_{AC}$ turns out to be Gaussian for each pair
$(s,s^{\prime})$, whereas the full reduced density matrix is not. When
specialised to the ground state of critical chains, these notations make it
manifest the correspondence pointed out in Ref. Coser2016Spin , according to
which
$(-1)^{4\vec{\varepsilon}\cdot\vec{\delta}}\mathcal{F}_{f}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](A,B,C)$ is the partition function of the
underlying fermionic quantum field theory on the Riemann surface associated
with the Rényi-$\alpha$ entropy, with given boundary conditions along the
cycles of the canonical homology basis. We now discuss this identification.

#### Spin structures: quantum field theory.

In QFT, the moments $\mathrm{tr}\rho_{A\cup C}^{\alpha}$ of the reduced
density matrix of disjoint blocks can be interpreted in terms of partition
functions on Riemann surfaces $\mathcal{R}_{\alpha}$ with genus $g=\alpha-1$.
We refer the interested reader to Refs Calabrese2009Entanglement ;
Calabrese2011Entanglement ; Headrick2013Bose ; Coser2014OnRenyi ;
Coser2016Spin ; Calabrese2011Entanglement  for more detailed presentations of
such surfaces; here we mention just their main properties. The Riemann surface
$\mathcal{R}_{\alpha}$ is given by the complex algebraic curve
$(y,z)\in\mathbb{C}^{2}$, with
$y^{\alpha}=(z-u_{1})(z-u_{2})[(z-v_{1})(z-v_{2})]^{\alpha-1}$. A key object
characterising a Riemann surface is the so-called period matrix. Since it
depends on the homology basis, we start by introducing such basis. That is a
set $\\{a_{j},b_{j};j=1,2,\ldots,g\\}$ of closed oriented curves that cannot
be contracted to a point, with the following intersection properties. Defining
the intersection number $h\circ\tilde{h}$ of two oriented curves $h,\tilde{h}$
as their number of intersections, with the orientation taken into account, the
cycles satisfy $a_{i}\circ a_{j}=b_{i}\circ b_{j}=0$ and $a_{i}\circ
b_{j}=-b_{i}\circ a_{j}=\delta_{ij}$. In this work we use the same canonical
cohomology basis as Ref. Coser2016Spin  (and also Refs Alvarez-Gaume1986 ;
Verlinde1987 ), shown, e.g., in Fig. 6 for the case $\alpha=4$. More details
can be found in Refs Coser2016Spin ; Coser2014OnRenyi . The period matrix for
this basis is given in Eq. (2.5). We remark that it is purely imaginary, a
property that is lost with more than two disjoint blocks Headrick2013Bose .
For $\alpha=2$ the Riemann surface is a torus and the period matrix reduces to
the so-called modular parameter $\tau=\Omega_{11}=i\beta_{1/2}(x)$, where
$\beta$ is defined in (2.5). The modular parameter $\tau$ can also be
expressed in terms of the complete elliptic integral $K$, as reported in
Figure 1.

For fermionic models on higher genus Riemann surfaces, the boundary conditions
along the cycles of the canonical homology basis can be either periodic
(Ramond) or antiperiodic (Neveu-Schwarz). The set of boundary conditions over
all cycles provides the spin structure of the model on the underlying Riemann
surface. We use standard notations according to which $\vec{\varepsilon}$
denotes the vector of the boundary conditions along $a$-cycles
($\varepsilon_{i}=0$ for antiperiodic and $\varepsilon_{i}=1/2$ for periodic
boundary conditions), whereas $\vec{\delta}$ is the vector of the boundary
conditions along $b$-cycles ($\delta_{i}=0$ for antiperiodic and
$\delta_{i}=1/2$ for periodic boundary conditions). It is not accidental that
the same symbols $\vec{\varepsilon}$ and $\vec{\delta}$ appear in (3.18). Such
a connection was established in Ref. Coser2016Spin , which showed that the
elements of the pseudo-RDMs introduced in (3.17) can be expressed in the
fermionic coherent state path integral formalism based on Grassmann variables
Kleinert2009 . Specifically, two Grassmann variables
$\zeta_{\ell},\zeta_{\ell}^{\dagger}$ are associated with each lattice site
$\ell$ and are then promoted to fields $\zeta(x),\zeta^{\dagger}(x)$ in the
continuum limit. Let us assume, for the sake of simplicity, that the system is
described by a single fermionic field $\bm{\mathbf{\psi}}(x)$, satisfying the
standard anticommutation relations
$\\{\bm{\mathbf{\psi}}(x),\bm{\mathbf{\psi}}^{\dagger}(y)\\}=\bm{\mathbf{I}}\delta(x-y)$.
The coherent states $\ket{\zeta}$ then satisfy
$\bm{\mathbf{\psi}}(x)\ket{\zeta}=\zeta(x)\ket{\zeta}$,
$\bra{\zeta}\bm{\mathbf{\psi}}^{\dagger}(x)=\bra{\zeta}\zeta^{\dagger}(x)$. We
also use the notation
$\ket{\zeta}=\ket{\zeta_{A},\zeta_{B},\zeta_{C},\zeta_{\overline{ABC}}}$ to
specify different subsystems separately, i.e. $\zeta(x)=\zeta_{A}(x)$ for
$x\in A$ etc. Then the elements of the pseudo RDMs can be written as

$\bra{\zeta_{A},\zeta_{C}}\rho_{AC}^{s,s^{\prime}}\ket{\eta_{A},\eta_{C}}=e^{-\zeta_{A}^{\dagger}\eta_{A}-\zeta_{C}^{\dagger}\eta_{C}}\\\
\int\mathcal{D}\zeta_{B}^{\dagger}\mathcal{D}\zeta_{B}\mathcal{D}\zeta_{\overline{ABC}}^{\dagger}\mathcal{D}\zeta_{\overline{ABC}}\braket{\zeta_{A},-s\zeta_{C},-s^{\prime}\zeta_{B},-\zeta_{\overline{ABC}}}{0}\braket{0}{\eta_{A},-s\eta_{C},\zeta_{B},\zeta_{\overline{ABC}}}\,,$
(3.21)

where $\zeta^{\dagger}_{A}\eta_{A}=\sum_{\ell\in
A}\zeta_{\ell}^{\dagger}\eta_{\ell}$.

In the continuum limit the overlaps $\braket{\zeta}{0}$ and
$\braket{0}{\zeta^{\prime}}$ can be represented as path integrals in the lower
and upper half plane respectively. Ref. Coser2016Spin  has proceeded along
these lines and reinterpreted Eq. (3.19) in terms of partition functions.
Namely,

$\mathrm{tr}\left[\prod\limits_{k=1}^{\alpha-1}\rho_{AC}^{s_{k}[\vec{\varepsilon}],s_{k}^{\prime}[\vec{\delta}]}\right]=Z_{\vec{\varepsilon},\vec{\delta}}=\int_{\vec{\varepsilon},\vec{\delta}}\mathcal{D}^{\dagger}\psi\mathcal{D}\psi\
e^{-S_{E}[\psi^{\dagger},\psi]},$ (3.22)

where $S_{E}=\int\mathcal{L}[\psi^{\dagger},\psi]dxd\tau$ is the Euclidean
action of the fermionic QFT description of the model, $\psi(x,\tau)$,
$\psi^{\dagger}(x,\tau)$ are Grassmann fields and the path integral is
performed on the Riemann surface $\mathcal{R}_{\alpha}$, with spin structure
$\vec{\varepsilon}$, $\vec{\delta}$ for the previously introduced canonical
homology basis. When there are multiple fermionic fields, such as the ones
constituting the Dirac spinor, the previous considerations generalize
straightforwardly and one works then with multiple fields with correlated spin
structures. Ref. Coser2016Spin  has in this way identified the spin structures
of the free Dirac and Majorana QFTs with the joint moments of the pseudo-RDMs
defined in Ref. Fagotti2010disjoint . We remark that, while both Ref.
Coser2016Spin  and Ref. Fagotti2010disjoint  focused on non-interacting spin
chains, the presented arguments are in fact general. In this work, in turn, we
will apply the same picture in the presence of interactions.

$a_{1}$$b_{1}$$a_{2}$$b_{2}$$a_{3}$$b_{3}$ Figure 6: The Riemann surface
arising in the computation of the Rényi-$4$ entropy of two disjoint blocks is
represented through the cyclic joining of four sheets. The lower edge of a
slit (black) should be identified with the upper edge of a slit above
(magenta), in a cyclic way. The canonical homology basis used in this work is
indicated. It is the same one as in Ref. Coser2016Spin .

###  3.2 Conformal critical systems

In this section we briefly review the behavior of the entanglement entropies
in critical spin chains described by a conformal field theory, and the
corresponding formalism. We focus on ground states in the thermodynamic limit.
More details about the formalism for computing the entanglement entropies in
QFT can be found in reviews Calabrese2009Entanglement ; Rangamani2017 .

When the subsystem is a single interval $A=(u,v)$, the moments of the reduced
density matrix for a criticial system described by a CFT with central charge
$c$ are given by Holzhey1994Geometric ; Calabrese2004Entanglement

$\mathrm{tr}\rho_{A}^{\alpha}=c_{\alpha}\left(\frac{v-u}{a}\right)^{-\frac{c}{6}\left(\alpha-\frac{1}{\alpha}\right)},$
(3.23)

where $a$ is the lattice spacing and $c_{\alpha}$ is a non-universal constant.
When the subsytem consists of two disjoint intervals, $A=(u_{1},v_{1})$ and
$C=(u_{2},v_{2})$, Refs Calabrese2009Entanglement1 ; Furukawa2009Mutual ;
Caraglio2008Entanglement  found

$\mathrm{tr}\rho_{A\cup
C}^{\alpha}=c_{\alpha}^{2}\left[\frac{a^{2}(u_{2}-u_{1})(v_{2}-v_{1})}{(v_{1}-u_{1})(v_{2}-u_{1})(v_{2}-u_{2})(u_{2}-v_{1})}\right]^{\frac{c}{6}\left(\alpha-\frac{1}{\alpha}\right)}\mathcal{F}(x),$
(3.24)

where $c_{\alpha}$ is the same non-universal constant as before, $x$ is the
cross ratio (1.4), which can be written as

$x=\frac{(v_{1}-u_{1})(v_{2}-u_{2})}{(u_{2}-u_{1})(v_{2}-v_{1})}\;,$ (3.25)

and $\mathcal{F}(x)$ is a universal function of the cross ratio that is model
dependent, in that it depends on the full operator content of the CFT. A very
convenient way of expressing the result for disjoint blocks is through the
tripartite information (1.3) of subsystems $A,B,C$, where $B$ is the interval
between $A$ and $C$,

$I_{3}^{(\alpha)}(A,B,C)=I_{3}^{(\alpha)}(x)=\frac{1}{\alpha-1}\log\mathcal{F}_{\alpha}(x).$
(3.26)

Thus, the universal part of the entanglement entropy of disjoint blocks is
determined by the tripartite information of three adjacent blocks.

The relation between the tripartite information of $A,B,C$ and the mutual
information $I_{2}^{(\alpha)}(A,C)=S_{\alpha}(A)+S_{\alpha}(C)-S_{\alpha}(AC)$
follows from Eqs (1.2) and (1.3)

$I_{3}^{(\alpha)}(A,B,C)=\frac{c}{6}\left(1+\frac{1}{\alpha}\right)\log(1-x)+I_{2}^{(\alpha)}(A,C)\;.$
(3.27)

In the limit $x\to 0$, corresponding to a large separation $|B|$ with respect
to the lengths of $A$ and $C$, the mutual information is expected to vanish.
Since the first term on the right hand side vanishes as well, we obtain
$I_{3}^{(\alpha)}(x)=0$ for $x\to 0$. This accounts for the choice of
normalization $\mathcal{F}_{\alpha}(0)=1$.

We mention that the result (3.23) for the single block can be interpreted as a
two-point function
$\mathrm{tr}\rho_{A}^{\alpha}=\braket{\mathcal{T}_{\alpha}(u)\overline{\mathcal{T}}_{\alpha}(v)}$
of branch-point twist fields Calabrese2004Entanglement ; Cardy2008Form  acting
at $u$ and $v$. The twist fields behave as spinless primary operators with
scaling dimension $\Delta_{\alpha}=\tfrac{c}{12}(\alpha-\tfrac{1}{\alpha})$.
Similarly, the result (3.24) can be interpreted as a four-point function
$\mathrm{tr}\rho_{A\cup
C}^{\alpha}=\braket{\mathcal{T}_{\alpha}(u_{1})\overline{\mathcal{T}}_{\alpha}(v_{1})\mathcal{T}_{\alpha}(u_{2})\overline{\mathcal{T}}_{\alpha}(v_{2})}$
and its form is dictated by global conformal invariance.

A general expectation is that the tripartite information $I_{3}^{(\alpha)}(x)$
(and the function $\mathcal{F}(x)$) exhibits the crossing symmetry
$x\leftrightarrow 1-x$. One way to see this is to consider a finite periodic
system of length $L$. There, the lengths $\ell$ entering the cross ratio are
replaced by chord lengths $\tfrac{L}{\pi}\sin(\tfrac{\ell\pi}{L})$, so that
the cross ratio becomes

$x=\frac{\sin\left(\pi\frac{v_{1}-u_{1}}{L}\right)\sin\left(\pi\frac{v_{2}-u_{2}}{L}\right)}{\sin\left(\pi\frac{u_{2}-u_{1}}{L}\right)\sin\left(\pi\frac{v_{2}-v_{1}}{L}\right)}\;.$
(3.28)

One has then

$I_{3}^{(\alpha)}(x)=I_{3}^{(\alpha)}(A,B,C)=I_{3}^{(\alpha)}(B,C,\overline{ABC})=I_{3}^{(\alpha)}(1-x)\,,$
(3.29)

where the second equality is a consequence of the property that in pure states
the entropy of a subsystem $X$ is equal to the entropy of its complement
$\bar{X}$, and the third equality stems from the property that the cross ratio
for disjoint intervals $B$ and $\overline{ABC}$ is equal to $1-x$.

###  3.3 Noncritical systems

#### Disordered phases.

In a state with clustering properties and a finite correlation length, at the
leading order in the lengths, the entanglement entropy of a subsystem of spins
depends only on the number of subsystem’s boundaries. This is known as “area
law” Hastings2007  and is sufficient to conclude that the tripartite
information is asymptotically zero. One could still wonder whether, in the
presence of spin-flip symmetry, something nontrivial could be found by
singling out the contribution coming from a given spin structure in the
underlying fermionic theory. We argue that the answer is no. The spin
structure in that case was briefly addressed in Ref. Fagotti2010disjoint
through a numerical analysis in the paramagnetic phase of the quantum XY
model. It was inferred

$\mathcal{F}_{f}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](x)=1\,.$ (3.30)

In Appendix A we lift this numerical observation to a prediction for any state
satisfying the aforementioned properties. It follows that the tripartite
information vanishes, $I_{3}^{(\alpha)}=0$.

#### Ordered phases.

Our proof of (3.30) relies on clustering and spin-flip symmetry; the only way
to circumvent (3.30) in a gapped system is therefore by breaking one of the
two. As we discussed before, spin flip is however important, if not necessary,
in order to give sense to the underlying fermionic subsystems, therefore we do
not allow ourselves to give it up in the representation through Jordan-Wigner
fermions. We thus examine symmetric ground states that break clustering. We
assume that the state is a cat state of the form

$\ket{\Psi}=\frac{1+(-1)^{\bm{\mathbf{F}}}}{\sqrt{2}}\ket{\psi}\,,$ (3.31)

where $\ket{\psi}$ has clustering properties and is macroscopically different
from $(-1)^{\bm{\mathbf{F}}}\ket{\psi}$. Ref. Fagotti2010disjoint  considered
also this case by numerically studying the ferromagnetic phase of the XY
model; specifically, it was inferred

$\mathcal{F}_{f}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](x)=(-1)^{4\vec{\varepsilon}\cdot\vec{\delta}}$
(3.32)

In Appendix A we provide arguments that this result should be expected in
general. It follows that the tripartite information in this case is
$I_{3}^{(\alpha)}=\log 2$.

#### Scaling close to discontinuous phase transitions.

In the rest of the section we investigate a limit in which the second class of
states approaches the former. Having in mind the physical interpretation of
the two classes as ground states of systems in disordered and ordered phases,
one could naively expect that, preserving the symmetry (spin flip), such a
limit would correspond to a quantum phase transition in which the entanglement
entropies exhibit logarithmic singularities (the reader could, for example,
think of the effect of changing the magnetic field in the quantum Ising
model). This is however avoidable in a presence of a point with a
discontinuous phase transition Franchini2007Ellipses ; Ercolessi2011Essential
, opening the door to the existence of states that satisfy the area law but
that exhibit a nontrivial tripartite information. We consider, in particular,
the following cat state

$\ket{\bar{\theta}}=\tfrac{\ket{\theta}+\ket{-\theta}}{\sqrt{2}}\,,$ (3.33)

where

$\ket{\theta}=e^{i\frac{\theta}{2}\sum_{\ell}\sigma_{\ell}^{y}}\ket{\Uparrow}$

and $\Uparrow$ denotes a string of all spin up. A reader familiar with the
quantum XY model Lieb1961 ; Barouch1971 , defined by the Hamiltonian

$\bm{\mathbf{H}}_{\mathrm{XY}}=\sum_{\ell}\left[(1+\gamma)\bm{\mathbf{\sigma}}_{\ell}^{x}\bm{\mathbf{\sigma}}_{\ell+1}^{x}+(1-\gamma)\bm{\mathbf{\sigma}}_{\ell}^{y}\bm{\mathbf{\sigma}}_{\ell+1}^{y}+2h\bm{\mathbf{\sigma}}_{\ell}^{z}\right]\;,$
(3.34)

might have recognised this as the symmetric ground state in a special region
of the phase diagram within its (anti)ferromagnetic phase Kurmann1982 ;
Muller1985Implications ; Franchini2017 , with
$\cos^{2}(\theta)=(1-\gamma)/(1+\gamma)$. The limit $\theta\rightarrow 0$
brings the cat state into a symmetric product state, which can be interpreted
as a ground state in a disordered phase (it is, indeed, the ground state in
the paramagnetic phase of the isotropic XY model).

First we note that we are free to replace the coherent superposition by an
incoherent one: the system is in the thermodynamic limit and $\ket{\theta}$ is
macroscopically different from $\ket{-\theta}$ for every $\theta\neq 0$. The
spin representation of the fermionic pseudo-RDM, (3.17), reads

$\rho^{s,s^{\prime}}=(-\cos\theta)^{\frac{1-s}{2}|B|}\rho^{(\theta;+)}_{A}\rho^{(\theta;+)}_{C}+s^{\prime}(-\cos\theta)^{\frac{1+s}{2}|B|}\rho^{(\theta;-)}_{A}\rho^{(\theta;-)}_{C}$
(3.35)

where

$\rho_{A}^{(\theta;s)}=\frac{e^{i\frac{\theta}{2}\sum_{\ell\in
A}\sigma_{\ell}^{y}}\ket{\Uparrow_{A}}\bra{\Uparrow_{A}}e^{-i\frac{\theta}{2}\sum_{\ell\in
A}\sigma_{\ell}^{y}}+se^{-i\frac{\theta}{2}\sum_{\ell\in
A}\sigma_{\ell}^{y}}\ket{\Uparrow_{A}}\bra{\Uparrow_{A}}e^{i\frac{\theta}{2}\sum_{\ell\in
A}\sigma_{\ell}^{y}}}{2}$ (3.36)

and, in particular, $\rho_{A}^{(\theta;+)}$ is equivalent to the reduced
density matrix of a spin block. The spin structure can be easily computed, and
we find

$\mathcal{F}_{f}\left[\begin{smallmatrix}\varepsilon\\\
\delta\end{smallmatrix}\right](A,B,C)=\tfrac{(1+\cos^{2(|A|+|B|)}(\theta))(1+\cos^{2(|B|+|C|)}(\theta))}{1+\cos^{2(|A|+|B|+|C|)}(\theta)}\\\
\Bigl{[}\tfrac{1}{1+\cos^{2(-1)^{2\varepsilon}|B|}(\theta)}+\tfrac{(-1)^{2\delta}\tanh(|A|\log\cos\theta)\tanh(|C|\log\cos\theta)}{1+\cos^{-2(-1)^{2\varepsilon}|B|}(\theta)}\Bigr{]}$
(3.37)

Consistently with (3.32), for any nonzero $\theta$ this approaches
$(-1)^{4\epsilon\delta}$ in the limit of infinite lengths. If however we take
that limit while scaling also the angle so that $\theta=\frac{y}{\sqrt{|B|}}$
with $y$ finite, we get

$\mathcal{F}_{f}\left[\begin{smallmatrix}\varepsilon\\\
\delta\end{smallmatrix}\right](A,B,C)\xrightarrow{|A|\sim|B|\sim|C|\sim\theta^{-2}\rightarrow\infty}\tfrac{(1+e^{-(\chi_{1}+1)y^{2}})(1+e^{-(\chi_{2}+1)y^{2}})}{1+e^{-(\chi_{1}+\chi_{2}+1)y^{2}}}\\\
\Bigl{[}\tfrac{1}{1+e^{-(-1)^{2\varepsilon}y^{2}}}+\tfrac{(-1)^{2\delta}\tanh(\frac{\chi_{1}}{2}y^{2})\tanh(\frac{\chi_{2}}{2}y^{2})}{1+e^{(-1)^{2\varepsilon}y^{2}}}\Bigr{]}$
(3.38)

where $\chi_{1}=\frac{|A|}{|B|}$ and $\chi_{2}=\frac{|C|}{|B|}$. This means
that, for any arbitrarily large subsystems, we can find angles $\theta$ for
which the spin structure is different from the basic values (3.30) and (3.32).
This will be exploited in Section 6 to uncover the spin structure of the
Heisenberg XXZ model.

##  4 Heisenberg XXZ model

###  4.1 Quantum field theory description

The XXZ spin chain (2.1) is critical in a range of parameters, e.g., in the
absence of magnetic field $h$ the critical region is $\Delta\in[-1,1]$. Its
low energy and large distance properties admit an effective Luttinger liquid
description, developed in Refs Luther1975 ; Haldane1980 . For a pedagogical
introduction to the effective quantum field theory description in terms of
Thirring and free compactified boson Lagrangian we refer the reader to Ref.
Affleck1988LesHouches .

The identification of the QFT is based on the Jordan-Wigner transformation. It
is first established perturbatively and is then lifted to a non-perturbative
level by fixing the QFT parameters phenomenologically, through the comparison
bewteen QFT predictions and exact lattice results. Specifically, the non-
interacting model is identified with a free Dirac fermion and the interacting
part, quartic in the Jordan-Wigner fermions, is treated as a small
perturbation. This yields the massless Thirring model, whose Lagrangian in the
Euclidean space reads

$\mathcal{L}=\bar{\psi}\gamma^{\mu}\partial_{\mu}\psi-\frac{1}{2}\pi\lambda\
\bar{\psi}\gamma^{\mu}\psi\bar{\psi}\gamma_{\mu}\psi,$ (4.1)

where the Dirac spinors have two components and the Dirac matrices satisfy
$\\{\gamma_{\mu},\gamma_{\nu}\\}=\delta_{\mu,\nu}$ (e.g.
$\gamma_{0}=\sigma^{y},\ \gamma_{1}=\sigma^{x}$). Bosonization is then used to
map the Thirring model to the free boson compactified on a circle, given by
the Lagrangian

$\mathcal{L}=\frac{R^{2}}{8\pi}\partial_{\mu}\varphi\partial^{\mu}\varphi,$
(4.2)

where the bosonic fields are periodic in the target space with period $2\pi$,
i.e. $\varphi\equiv\varphi+2\pi$, and the compactification radius $R$ is
related to the Thirring coupling constant through $R^{2}=1+\lambda$. Note that
throughout the literature different conventions for the compactification
radius are used. In the one we use the T-duality reads $R\leftrightarrow 2/R$.

Next, symmetry arguments are used to conclude that there are no other relevant
Lorentz-invariant interactions in a range of $R$, while the various non-
renormalizible interactions can only lead to a renormalization of $R$. The
compactification radius is finally determined phenomenologically, by comparing
the lattice result for some quantity, typically obtained within Bethe ansatz,
with the field theoretical one. Various quantities can be used to this
purpose, such as some critical exponents Luther1975 , finite-size spectrum
Haldane1980 ; Alcaraz1987 ; Lukyanov1998 , compressibility Giamarchi2003book
and susceptibility Eggert1994 . For zero magnetic field $h=0$ the
identification is

$\eta=\frac{1}{2}R^{2}=\frac{1}{2}(1+\lambda)=\frac{1}{2K_{L}}=\frac{\arccos(-\Delta)}{\pi}\;,$
(4.3)

where we have also indicated the relation with the Luttinger liquid parameter
$K_{L}$ (see e.g. Giamarchi2003book ). The exponent $\eta$ appears in the
spin-correlation functions, whose large distance behavior for $h=0$ is given
by Luther1975

$\braket{\bm{\mathbf{\sigma}}^{x}_{0}\bm{\mathbf{\sigma}}^{x}_{\ell}}\sim(-1)^{\ell}\frac{1}{|\ell|^{\eta}},\qquad\braket{\bm{\mathbf{\sigma}}^{z}_{0}\bm{\mathbf{\sigma}}^{z}_{\ell}}\sim\frac{1}{\ell^{2}}+(-1)^{\ell}\frac{1}{|\ell|^{1/\eta}}\;.$
(4.4)

For a non-zero magnetic field the last equality in Eq. (4.3) does not hold,
nor Eq. (4.4) for the spin correlation functions. In that case $\eta$ is
obtained numerically by solving integral equations Cabra1998  or through
perturbative methods Granet2019 .

###  4.2 The free boson compactified on a circle

The Rényi-$\alpha$ tripartite information for the free compactified boson was
first computed for $\alpha=2$ in Ref. Furukawa2009Mutual  and then was
generalized to integer $\alpha\geq 2$ in Ref. Calabrese2009Entanglement1 . The
result reads

$I_{3}^{(\alpha)}(x)=\frac{1}{\alpha-1}\log\frac{\Theta(\vec{0}|\eta\Omega)\Theta(\vec{0}|\Omega/\eta)}{\left[\Theta(\vec{0}|\Omega)\right]^{2}},$
(4.5)

where $\Omega$ is the period matrix Calabrese2009Entanglement1  of the Riemann
surface with genus $\alpha-1$, given in Eq. (2.5), and
$\Theta\equiv\Theta[\vec{0}]$ is the Riemann-Siegel theta function defined in
Eq. (2.6). The result (4.5) was derived in Ref. Calabrese2009Entanglement1
using older results about four-point correlation functions of twist fields
Dixon1987 . The result was rederived in Ref. Coser2014OnRenyi , where also the
case of more than two intervals is discussed, using older results on partition
functions on higher genus Riemann surfaces Zamolodchikov1987 ; Alvarez-
Gaume1986 ; Alvarez-Gaume1987 ; Verlinde1987 ; Dijkgraaf1988 . The formula
(4.5) is manifestly invariant under $\eta\to 1/\eta$, i.e. $R\to 2/R$.
Although not apparent from the expression above, the result also respects the
crossing symmetry $x\leftrightarrow 1-x$.

The analytical continuation of (4.5) and, accordingly, the limit
$\alpha\rightarrow 1^{+}$ are not known in general. However, it is possible to
perform the limit $\alpha\rightarrow 1^{+}$ for small $x$. The small $x$
expansion for $\eta\neq 1$ reads Calabrese2011Entanglement

$I_{3}^{(\alpha)}(x)=\frac{s\left(\alpha;\mathcal{A}_{\textrm{CB}}\right)}{\alpha-1}\left(\frac{x}{4\alpha^{2}}\right)^{\mathcal{A}_{\textrm{CB}}}+o\left(x^{\mathcal{A}_{\textrm{CB}}}\right),\qquad\eta\neq
1,\quad\textrm{as }x\to 0,$ (4.6)

where

$\mathcal{A}_{\textrm{CB}}=\mathrm{min}[\eta,1/\eta]$ (4.7)

and we have defined the function

$s(\alpha;\mathcal{A})=\alpha\sum_{j=1}^{\alpha-1}\left[\sin\left(\pi\frac{j}{\alpha}\right)\right]^{-2\mathcal{A}}\;,\qquad\alpha=2,3,4,\ldots$
(4.8)

Ref. Calabrese2011Entanglement  also obtained the derivative at $\alpha=1$ of
the analytical continuation of $s(\alpha;\mathcal{A})$ to complex $\alpha$,

$s_{\mathrm{vN}}(\mathcal{A})\equiv\frac{d}{d\alpha}s(\alpha;\mathcal{A})\Big{|}_{\alpha=1}=\frac{\sqrt{\pi}\Gamma(\mathcal{A}+1)}{2\Gamma(\mathcal{A}+\frac{3}{2})},$
(4.9)

given also in Eq. (2.8). This result can be used to obtain the small $x$
expansion of the von Neumann tripartite information,

$I_{3}^{(\mathrm{vN})}(x)=\frac{\sqrt{\pi}\Gamma(\mathcal{A}_{\textrm{CB}}+1)}{2\Gamma(\mathcal{A}_{\textrm{CB}}+\frac{3}{2})}\left(\frac{x}{4}\right)^{\mathcal{A}_{\textrm{CB}}}+o\left(x^{\mathcal{A}_{\textrm{CB}}}\right),\quad\textrm{as
}x\to 0$ (4.10)

At the self-dual point $\eta=1$ the tripartite information (4.5) is exactly
zero for any $\alpha$, and therefore also for $\alpha=1$.

The formula (4.5) has been checked numerically in the XXZ chain for the lowest
values of $\alpha$ using exact diagonalization Furukawa2009Mutual  and tree
tensor networks Alba2011Entanglement . We stress that it would not be correct
to identify $\eta$ with $1/\eta$ in the chain. Namely, although the tripartite
information in Eq. (4.5) is invariant under $\eta$ and $1/\eta$, some
quantities, like the spin correlation functions in Eq. (4.4), are not. Note
that the identification (4.3) for the XXZ chain without magnetic field makes
it evident that $\eta\in(0,1]$ so the duality $\eta\leftrightarrow 1/\eta$
cannot be addressed there.

###  4.3 Thirring model on a Riemann surface

Here we point out that the tripartite information in the XXZ chain can also be
computed within the fermionic approach, by computing (3.19) for each
$\vec{\varepsilon},\vec{\delta}$ separately and then summing them. As
discussed in section 3.1, the last term in (3.19) in the fermionic field
theory becomes a partition function on the higher genus Riemann surface with
the specified spin structure. Although the arguments are general, the method
has been applied so far only in non-interacting spin chains.

To derive the tripartite information in the XXZ chain in the fermionic picture
we need to evaluate the partition functions of the massless Thirring model on
a higher genus closed Riemann surface for different spin structures. Several
papers Freedman1988Thirring ; Freedman1989Thirring ; Wu1989Determinants ;
Sachs1996  have already tackled with this problem, by manipulating fermionic
path integrals, and the results display inconsistencies about factors and the
sign of the partition functions, as e.g. commented in Ref. Sachs1996 . In this
work we fix such ambiguities by requiring that our results reproduce the known
predictions for the entanglement entropies of disjoint blocks in the bosonic
picture Calabrese2009Entanglement1 . In this way we fix some sign
discrepancies but we do not investigate overall factors common to all spin
structures, which are irrelevant to our purposes.

We conclude that, for purely imaginary period matrix $\Omega$, we have

$Z(\vec{\varepsilon},\vec{\delta};\lambda)=(-1)^{4\vec{\varepsilon\cdot\delta}}Z(\vec{0},\vec{0};0)\
f_{\alpha}(\lambda)\ \frac{\Theta\left[\begin{matrix}\vec{E}\\\
\vec{F}\end{matrix}\right]\left(\vec{0}|\tilde{\Omega}\right)}{\left[\Theta(\vec{0}|\Omega)\right]^{2}}\;,$
(4.11)

where $Z(\vec{0},\vec{0};0)$ is the partition function for the free Dirac
fermion with all antiperiodic boundary coniditions, $f_{\alpha}$ are unknown
functions of $\lambda$ (neglected in Fig. 1) such that $f_{\alpha}(0)=0$,
$\tilde{\Omega}$ is a $2(\alpha-1)\times 2(\alpha-1)$ matrix defined by

$\tilde{\Omega}=\Bigl{[}\bigl{(}\tfrac{1}{4\eta}+\eta\bigr{)}\mathrm{I}+\bigl{(}\tfrac{1}{4\eta}-\eta\bigr{)}\sigma^{x}\Bigr{]}\otimes\Omega\;,\quad\eta=\frac{1}{2}(1+\lambda),$
(4.12)

and the vectors

$\vec{E}=\begin{pmatrix}1\\\
-1\end{pmatrix}\otimes\vec{\varepsilon}\,,\qquad\vec{F}=\begin{pmatrix}1\\\
1\end{pmatrix}\otimes\vec{\delta}$ (4.13)

carry the information about the spin structure. The formula (4.11) differs
from the result of Refs Freedman1988Thirring ; Freedman1989Thirring ;
Wu1989Determinants  in the sign $(-1)^{4\varepsilon\cdot\delta}$. In the free
case $\lambda=0$, studied in Ref. Alvarez-Gaume1986 , the partition function
of odd spin structures vanishes so that sign is irrelevant and formula (4.11)
captures the relation between different spin structures

$\frac{Z(\vec{\varepsilon},\vec{\delta};0)}{Z(\vec{0},\vec{0};0)}=\Bigg{|}\Theta\left[\begin{matrix}\vec{\varepsilon}\\\
\vec{\delta}\end{matrix}\right]\left(\vec{0}|\Omega\right)/\Theta\left(\vec{0}|\Omega\right)\Bigg{|}^{2}\,.$
(4.14)

###  4.4 Spin structures

In this section we compute the function $\mathcal{F}_{f}$, defined in (3.19),
for each spin structure in the XXZ chain. We know that the partition function
$Z(\vec{0},\vec{0};0)$ is equal, up to a non-universal multiplicative
constant, to the trace in Eq. (3.24) with $c=1$ and $\mathcal{F}(x)=1$, where
the latter is a result of Ref. Casini2005 ; Casini2009reduced_density ;
Casini2009Entanglement . With this insight, and using also (3.23), we conclude

$\mathcal{F}_{f}\left[\begin{matrix}\vec{\varepsilon}\\\
\vec{\delta}\end{matrix}\right]=(-1)^{4\vec{\varepsilon}\cdot\vec{\delta}}\mathcal{C}_{\alpha}\frac{\Theta\left[\begin{matrix}\vec{E}\\\
\vec{F}\end{matrix}\right]\left(\vec{0}|\tilde{\Omega}\right)}{\left[\Theta(\vec{0}|\Omega)\right]^{2}}\;,$
(4.15)

for some non-universal constant $\mathcal{C_{\alpha}}$, that is independent of
the spin structure. The numerator in Eq. (4.15) is expressed in terms of a
$2(\alpha-1)\times 2(\alpha-1)$ matrix $\tilde{\Omega}$. It is possible to
express it in terms of $(\alpha-1)\times(\alpha-1)$ matrices, by resorting to
the identity

$\Theta\left[\begin{matrix}\vec{E}\\\
\vec{F}\end{matrix}\right]\left(\vec{0}|\tilde{\Omega}\right)=\sum_{\vec{\mu}\in\\{0,\frac{1}{2}\\}^{\alpha-1}}(-1)^{4\vec{\mu}\cdot\vec{\delta}}\Theta\left[\begin{matrix}\vec{\mu}\\\
\vec{0}\end{matrix}\right]\left(\vec{0}|\tfrac{1}{\eta}\Omega\right)\Theta\left[\begin{matrix}\vec{\varepsilon}+\vec{\mu}\\\
\vec{0}\end{matrix}\right]\left(\vec{0}|4\eta\Omega\right)\,,$ (4.16)

proven in Appendix B.

The tripartite information (3.18) is obtained by summing over all spin
structures. The constant $\mathcal{C}_{\alpha}$ can be fixed by requiring that
the tripartite information vanishes at $x\to 0$, which is a consequence of
clustering (an equivalent trick was used in Ref. Calabrese2011Entanglement
for the Ising CFT). This fixes the constant to $\mathcal{C}_{\alpha}=1$. In
this way we obtain the result in Eq. (2.10). In Appendix B we show that
formula (B.1) reproduces the known result (4.5) obtained from studying the
free compactified boson.

#### Small $x$ expansion

The small $x$ expansion of (2.10) can be derived using the methods developed
in Ref. Calabrese2011Quantum  for the free compactified boson and Ising CFT.
Details of the procedure are presented in Appendix F. We are going to look
separately at three different classes of spin structures.

  1. 1)

Spin structures with $\vec{\varepsilon}=\vec{0}$
These are the only spin structures resulting in a non-zero value in the limit
$x\to 0$. The small $x$ expansion is

$\begin{split}&\mathcal{F}_{f}\left[\begin{matrix}\vec{0}\\\
\vec{\delta}\end{matrix}\right]=1+s(\alpha;1/\eta)\left(\frac{x}{4\alpha^{2}}\right)^{\tfrac{1}{\eta}}+s(\alpha;4\eta)\left(\frac{x}{4\alpha^{2}}\right)^{4\eta}-s(\alpha;1)\frac{x}{2\alpha^{2}}\\\
&+\sum\limits_{0\leq j_{1}<j_{2}\leq\alpha-1}\frac{4\
(-1)^{4\sum_{k=j_{1}+1}^{j_{2}}\delta_{k}}}{\left[\sin\left(\pi\frac{j_{2}-j_{1}}{\alpha}\right)\right]^{2\big{(}\tfrac{1}{4\eta}+\eta\big{)}}}\left(\frac{x}{4\alpha^{2}}\right)^{\tfrac{1}{4\eta}+\eta}+o\big{(}x^{\min(\tfrac{1}{\eta},4\eta,1,\tfrac{1}{4\eta}+\eta)}\big{)}\end{split}$
(4.17)

where $s$ is defined in Eq. (4.8).

  2. 2)

Spin structures with two “domain walls” in $\vec{\varepsilon}$ (imagining
$\varepsilon_{-1}=\varepsilon_{\alpha}=0$ at the boundaries), i.e.

$\vec{\varepsilon}^{\mathrm{T}}=\frac{1}{2}(\underbrace{0,0,\ldots,0}_{j_{1}},\underbrace{1,1,\ldots,1}_{j_{2}-j_{1}},0,0,\ldots,0)\;,\quad
0\leq j_{1}<j_{2}\leq\alpha-1\;.$ (4.18)

The small $x$ expansion is

$\begin{split}&\mathcal{F}_{f}\left[\begin{matrix}\vec{\varepsilon}\\\
\vec{\delta}\end{matrix}\right]\\\
&=\frac{2}{\left[\sin\left(\pi\frac{j_{2}-j_{1}}{\alpha}\right)\right]^{\tfrac{1}{2\eta}}}\left(\frac{x}{4\alpha^{2}}\right)^{\tfrac{1}{4\eta}}+(-1)^{4\vec{\varepsilon}\cdot\vec{\delta}}\frac{2}{\left[\sin\left(\pi\frac{j_{2}-j_{1}}{\alpha}\right)\right]^{2\eta}}\left(\frac{x}{4\alpha^{2}}\right)^{\eta}+o\left(x^{\min(\eta,\tfrac{1}{4\eta})}\right)\end{split}$
(4.19)

independently of $\vec{\delta}$.

  3. 3)

Spin structures with four or more “domain walls” in $\vec{\varepsilon}$, i.e.
all the remaining spin structures.
These spin structures exist only for $\alpha\geq 4$ and they are subleading
with respect to the spin-structures of type $1$,

$\mathcal{F}_{f}\left[\begin{matrix}\vec{\varepsilon}\\\
\vec{\delta}\end{matrix}\right]=O\left(x^{\min(2\eta,\tfrac{1}{2\eta})}\right)\;.$
(4.20)

The small $x$-expansion for the tripartite information in the XXZ chain, given
by Eq. (4.6), comes from the first term in (4.17) and the second term in
(4.19) (and taking the logarithm). Note that the second term in (4.19) is
subleading for $\eta>1/2$, but summing over spin structures kills the leading
term.

##  5 Fermionic XXZ and models with higher central charge

###  5.1 Fermionic XXZ chain

Figure 7: Ratio between the Rényi-$2$ tripartite information in the ground
state of the XXZ fermionic chain (2.2) and that in the XXZ spin chain (2.1) vs
$\Delta$ for $x=1/2$ and $h=0$. For $\Delta=0$, corresponding to non-
interacting chains, the tripartite information in the fermionic model is zero.
In its vicinity the tripartite information in the fermionic chain is much
smaller than the one in the corresponding spin chain.

In the previous section we have computed the contribution of each spin
structure to the tripartite information of the XXZ chain. From the result for
the spin structure $\vec{\varepsilon}=\vec{\delta}=\vec{0}$ (antiperiodic
boundary conditions on all cycles), we can readily compute the tripartite
information in the fermionic XXZ chain (2.2). In this way we obtain the result
presented in Eq. (2.4). In Fig. 7 the result is compared to the one for the
XXZ spin chain with $\alpha=2$. Prediction (2.4) is manifestly invariant under
$\eta\leftrightarrow 1/(4\eta)$, which is a symmetry around the free point
$\eta=1/2$. This symmetry should be contrasted to the symmetry
$\eta\leftrightarrow 1/\eta$ of the free compactified boson in Eq. (4.5). At
the free point $\eta=1/2$ we get $I_{3,\mathrm{f}}^{(\alpha)}(x)=0$, because
of the identity Headrick2013Bose

$|\Theta(\vec{0}|\Omega)|^{2}=\sum_{\vec{\varepsilon}\in\\{0,\frac{1}{2}\\}^{\alpha-1}}\Bigl{|}\Theta\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{0}\end{smallmatrix}\right]\left(\vec{0}|2\Omega\right)\Bigr{|}^{2}\,.$
(5.1)

The vanishing of the tripartite information in the non-interacting case, i.e.
in the Kitaev chain at its phase transition, is in agreement with the known
result for the free Dirac fermion Casini2005 ; Casini2009Entanglement ;
Casini2009reduced_density .

#### Small $x$ expansion.

Putting $\vec{\delta}=\vec{0}$ in (4.17) and taking the logarithm, we get the
small $x$ expansion for the tripartite information in the fermionic XXZ chain,

$\begin{split}I_{3,\mathrm{f}}^{(\alpha)}(x)=&\frac{1}{\alpha-1}\bigg{[}-s(\alpha;1)\frac{x}{2\alpha^{2}}+s(\alpha;\tfrac{1}{\eta})\left(\frac{x}{4\alpha^{2}}\right)^{\tfrac{1}{\eta}}+s(\alpha;4\eta)\left(\frac{x}{4\alpha^{2}}\right)^{4\eta}\\\
&+s(\alpha;\tfrac{1}{4\eta}+\eta)\left(\frac{x}{4\alpha^{2}}\right)^{\tfrac{1}{4\eta}+\eta}+o\big{(}x^{\min(1,\tfrac{1}{\eta},4\eta)}\big{)}\bigg{]}\end{split}$
(5.2)

Note that the term with power $\tfrac{1}{4\eta}+\eta$ contributes only for
$\eta=1/2$. It is instructive to single out the $\eta$-dependent leading term,

$I_{3,\mathrm{f}}^{(\alpha)}(x)=\begin{cases}s(\alpha;4\eta)\left(\frac{x}{4\alpha^{2}}\right)^{4\eta}+o(x^{4\eta}),&\eta<\tfrac{1}{4}\\\
-s(\alpha;1)\frac{x}{4\alpha^{2}}+o(x),&\eta=\tfrac{1}{4}\\\
-s(\alpha;1)\frac{x}{2\alpha^{2}}+o(x),&\tfrac{1}{4}<\eta<1,\eta\neq\tfrac{1}{2}\\\
-s(\alpha;1)\frac{x}{4\alpha^{2}}+o(x),&\eta=1\\\
s(\alpha;\tfrac{1}{\eta})\left(\frac{x}{4\alpha^{2}}\right)^{\tfrac{1}{\eta}}+o(x^{\tfrac{1}{\eta}}),&\eta>1\\\
\end{cases}$ (5.3)

It is now possible to compute the small $x$ expansion of the von Neumann
tripartite information using formula (4.9), derived in Ref.
Calabrese2011Entanglement  in the context of the spin chain. The result is
given in Eq. (2.7), or, singling out the leading term, by

$I_{3,\mathrm{f}}^{(\mathrm{vN})}(x)=\begin{cases}s_{\mathrm{vN}}(4\eta)\left(\frac{x}{4}\right)^{4\eta}+o(x^{4\eta}),&\eta<\tfrac{1}{4}\\\
-s_{\mathrm{vN}}(1)\frac{x}{4}+o(x),&\eta=\tfrac{1}{4}\\\
-s_{\mathrm{vN}}(1)\frac{x}{2}+o(x),&\tfrac{1}{4}<\eta<1,\eta\neq\tfrac{1}{2}\\\
-s_{\mathrm{vN}}(1)\frac{x}{4}+o(x),&\eta=1\\\
s_{\mathrm{vN}}(\tfrac{1}{\eta})\left(\frac{x}{4}\right)^{\tfrac{1}{\eta}}+o(x^{\tfrac{1}{\eta}}),&\eta>1\\\
\end{cases}$ (5.4)

Interestingly, while the tripartite information in the XXZ chain is positive
(or zero for $\eta=1$), we can see clearly that for small $x$ the tripartite
information in the fermionic XXZ chain is negative in the parameter region
$\eta\in[\tfrac{1}{4},\tfrac{1}{2})\cup(\tfrac{1}{2},1]$.

###  5.2 Jordan-Wigner deformations

There is a straightforward way to construct spin-chain models with higher
central charge starting from a massless system with local interactions. It is
sufficient to move operators apart in such a way that their distance is
multiplied by a given number $n$, e.g.,
$\sigma_{\ell}^{\alpha}\sigma_{\ell+1}^{\beta}\sigma_{\ell+3}^{\gamma}\rightarrow\sigma_{\ell}^{\alpha}\sigma_{\ell+n}^{\beta}\sigma_{\ell+3n}^{\gamma}$.
The resulting Hamiltonian is written as a sum of $\nu$ commuting extensive
operators acting on distinct sublattices and, appropriately rescaling the
lengths ($\ell\rightarrow\frac{\ell}{\nu}$), the low-energy physics is
described by the sum of independent CFTs. Such a theory is arguably
uninteresting.

If the Hamiltonian is invariant under a spin flip, a Jordan-Wigner (JW)
transformation allows one to reinterpret the spin-chain Hamiltonian as a
fermion one with local interactions. This opens the door to applying the same
procedure to the fermion Hamiltonian, which in the spin representation has the
effect of moving the spins apart, but additional spin operators in the
direction of the quantisation axis can appear. For example we have

$\sigma_{\ell}^{x}\sigma_{\ell+2}^{x}\xrightarrow{\nu}\sigma_{\ell}^{x}(\prod_{j=1\atop
j\neq\nu}^{2\nu-1}\sigma_{\ell+j}^{z})\sigma_{\ell+2\nu}^{x}\,,$

which are also known as “cluster interactions” Raussendorf2001 ; Doherty2009 ;
Son2011 ; Giampaolo2015 . Because of the spin structure, such a deformation
has now a nontrivial effect. This leads us to Hamiltonian (2.3), that we call
“Jordan-Wigner-$\nu$ deformations of the XXZ chain”. We will see in Sec. 6
that JW deformations characterise the spin structure almost completely.

Applying the JW transformation to the JW deformation of the XXZ chain, given
by Hamiltonian (2.3), we get

$\bm{\mathbf{H}}_{\mathrm{XXZ}}^{(\nu)}=2\sum_{\ell}\bm{\mathbf{c}}_{\ell+\nu}\bm{\mathbf{c}}^{\dagger}_{\ell}+\bm{\mathbf{c}}_{\ell}\bm{\mathbf{c}}^{\dagger}_{\ell+\nu}-2(\Delta+h)\bm{\mathbf{c}}^{\dagger}_{\ell}\bm{\mathbf{c}}_{\ell}+2\Delta\bm{\mathbf{c}}^{\dagger}_{\ell}\bm{\mathbf{c}}_{\ell}\bm{\mathbf{c}}^{\dagger}_{\ell+\nu}\bm{\mathbf{c}}_{\ell+\nu}$
(5.5)

The model is a sum of fermionic XXZ chain Hamiltonians (see Eq. (2.2)), on
sublattices $\nu\mathbb{Z}+j=\\{,\ldots,-\nu+j,j,\nu+j,2\nu+j\ldots\\}$, for
$j=0,1,\ldots,\nu-1$. The Hamiltonians are mutually commuting so they can be
diagonalized simultaneously. In particular, all the fermionic correlation
functions in the ground state factorize into those involving only fermions on
the same sublattice. This also implies the factorization of the function in
Eq. (4.15),

$\mathcal{F}_{f}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](A,B,C)=\prod_{j=1}^{\nu}\mathcal{F}_{f}^{(\mathrm{XXZ})}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](A|_{j},B|_{j},C|_{j}),$ (5.6)

where we use the notation $A|_{j}=A\cap(\nu\mathbb{Z}+j)$ for the intersection
between set $A$ and sublattice $j$. The factorization (5.6) is exact on the
lattice. In the limit of large subsystems we have a dependence on the
configuration only through the cross ratio, so

$\mathcal{F}_{f}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](x)=\left(\mathcal{F}_{f}^{(\mathrm{XXZ})}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](x)\right)^{\nu},$ (5.7)

Thus, the contribution of each spin structure to the tripartite information in
the Jordan-Wigner deformations is trivially related to the corresponding
contribution in the XXZ chain. By summing over spin structures, we get the
tripartite information in the Jordan-Wigner deformation of the XXZ chain,

$I_{3}^{(\alpha)}(x)=\frac{1}{\alpha-1}\log\left(\frac{1}{2^{\alpha-1}}\sum_{\vec{\varepsilon},\vec{\delta}\in\\{0,\frac{1}{2}\\}^{\alpha-1}}(-1)^{4\vec{\varepsilon}\cdot\vec{\delta}}\left(\mathcal{F}_{f}^{(\mathrm{XXZ})}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](x)\right)^{\nu}\right),$ (5.8)

where $\mathcal{F}_{f}^{(\mathrm{XXZ})}$ is given by Eq. (2.10). In the free
case $\eta=1/2$ the result reduces to

$I_{3}^{(\alpha;\nu)}(x)=\frac{1}{\alpha-1}\log\left(\frac{1}{2^{\alpha-1}}\sum_{\vec{\varepsilon},\vec{\delta}\in\\{0,\frac{1}{2}\\}^{\alpha-1}}(-1)^{4{\vec{\varepsilon}}\cdot{\vec{\delta}}}\left|\frac{\Theta\left[\begin{smallmatrix}\vec{\epsilon}\\\
\vec{\delta}\end{smallmatrix}\right](0|\Omega)}{\Theta(0|\Omega)}\right|^{2\nu}\right)\,,$
(5.9)

which could be deduced already from Refs Fagotti2012New ; Coser2016Spin .

#### Jordan-Wigner 2-deformation

Let us focus on the particular case $\nu=2$ and derive a simplified
expression. Writing explicitly the square as a product we have

$\begin{split}&\frac{1}{2^{\alpha-1}}\sum_{\vec{\varepsilon},\vec{\delta}\in\\{0,\frac{1}{2}\\}^{\alpha-1}}(-1)^{4\vec{\varepsilon}\cdot\vec{\delta}}\left(\mathcal{F}_{f}^{(\mathrm{XXZ})}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](x)\right)^{2}=\\\
&\frac{1}{2^{\alpha-1}}\sum_{\vec{\varepsilon},\vec{\delta}\in\\{0,\frac{1}{2}\\}^{\alpha-1}}(-1)^{4\vec{\varepsilon}\cdot\vec{\delta}}\sum_{\vec{\mu}_{1},\vec{\mu_{2}}\in\\{0,\frac{1}{2}\\}^{\alpha-1}}\prod_{j=1}^{2}(-1)^{4\vec{\mu}_{j}\cdot\vec{\delta}}\frac{\Theta\left[\begin{smallmatrix}\vec{\epsilon}+\vec{\mu}_{j}\\\
\vec{0}\end{smallmatrix}\right]\left(\vec{0}|\frac{1}{\eta}\Omega\right)\Theta\left[\begin{smallmatrix}\vec{\mu}_{j}\\\
\vec{0}\end{smallmatrix}\right]\left(\vec{0}|4\eta\Omega\right)}{\left[\Theta(\vec{0}|\Omega)\right]^{2}}\end{split}$
(5.10)

The sum over $\vec{\delta}$-dependent terms is

$\frac{1}{2^{\alpha-1}}\sum_{\vec{\delta}\in\\{0,\frac{1}{2}\\}^{\alpha-1}}(-1)^{4(\vec{\varepsilon}+\vec{\mu}_{1}+\vec{\mu}_{2})\cdot\vec{\delta}}=\begin{cases}1,&2(\vec{\mu}_{1}+\vec{\mu}_{2})=2\vec{\varepsilon}\quad(\mathrm{mod}\
2)\\\ 0,&\mathrm{otherwise}\end{cases}$ (5.11)

Since $\Theta[\vec{e}]=\Theta[\tfrac{1}{2}(2\vec{e}\ \mathrm{mod}\ 2)]$ for
$\vec{e}\in(\mathbb{Z}/2)^{2(\alpha-1)}$, after summing over $\vec{\delta}$ in
Eq. (5.10) we can replace $\vec{\varepsilon}$ by $\vec{\mu}_{1}+\vec{\mu}_{2}$
in the characteristic of the theta functions . For the same reason, vectors
$2\vec{\mu}_{2}+\vec{\mu}_{1}$ and $2\vec{\mu}_{1}+\vec{\mu}_{2}$ in the
characteristic can be replaced by $\vec{\mu}_{1}$ and $\vec{\mu}_{2}$
respectively. Grouping the factors that depend on $\vec{\mu}_{1}$ and those
that depend on $\vec{\mu}_{2}$ we recognize

$\frac{1}{2^{\alpha-1}}\sum_{\vec{\varepsilon},\vec{\delta}\in\\{0,\frac{1}{2}\\}^{\alpha-1}}(-1)^{4\vec{\varepsilon}\cdot\vec{\delta}}\left(\mathcal{F}_{f}^{(\mathrm{XXZ})}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](x)\right)^{2}=\left(\mathcal{F}_{f}^{(\mathrm{XXZ})}\left[\begin{smallmatrix}\vec{0}\\\
\vec{0}\end{smallmatrix}\right](x)\right)^{2}\;.$ (5.12)

Therefore, the tripartite information in the Jordan-Wigner 2-deformation is
the same as the fermionic one in the corresponding fermion model, which is
equal to twice the one in the fermionic XXZ chain, given by Eq. (2.4). This
result is expressed by Eq. (2.14).

#### Small $x$ expansion

For $\nu=2$ we have shown that the tripartite information is equal to twice
the one of the fermionic XXZ chain (see Eq. (2.14)). Therefore its small $x$
expansion can be read from the small $x$ expansion of the latter, given by Eq.
(5.2) and Eqs (5.3) for Rényi index $\alpha=2,3,\ldots$ and by Eqs (2.7) and
(5.4) for the von Neumann case.

For $\nu\geq 3$ the result can be obtained from the small $x$ expansion for
different spin structures in the XXZ chain (Eqs (4.17), (4.19) and (4.20)), by
raising them to power $\nu$, in accordance with Eq. (5.8). The result reads

$\displaystyle
I_{3}^{(\alpha;\nu=3)}(x)=\frac{1}{\alpha-1}\left[-3s(\alpha;1)\frac{x}{2\alpha^{2}}+3s(\alpha;\tfrac{1}{\eta})\left(\frac{x}{4\alpha^{2}}\right)^{\tfrac{1}{\eta}}+4s(\alpha;3\eta)\left(\frac{x}{4\alpha^{2}}\right)^{3\eta}+o\big{(}x^{\min(1,\tfrac{1}{\eta},3\eta)}\big{)}\right]\;,$
(5.13) $\displaystyle I_{3}^{(\alpha;\nu\geq
4)}(x)=\frac{\nu}{\alpha-1}\left[-s(\alpha;1)\frac{x}{2\alpha^{2}}+s(\alpha;\tfrac{1}{\eta})\left(\frac{x}{4\alpha^{2}}\right)^{\tfrac{1}{\eta}}+s(\alpha;4\eta)\left(\frac{x}{4\alpha^{2}}\right)^{4\eta}+o\big{(}x^{\min(1,\tfrac{1}{\eta},4\eta)}\big{)}\right]$
(5.14)

for the Rényi-$\alpha$ tripartite information, and

$\displaystyle
I_{3}^{(\mathrm{vN};\nu=3)}(x)=-3s_{\mathrm{vN}}(1)\frac{x}{2}+3s_{\mathrm{vN}}(\tfrac{1}{\eta})\left(\frac{x}{4}\right)^{\tfrac{1}{\eta}}+4s_{\mathrm{vN}}(3\eta)\left(\frac{x}{4}\right)^{3\eta}+o\big{(}x^{\min(1,\tfrac{1}{\eta},3\eta)}\big{)}\;,$
(5.15) $\displaystyle I_{3}^{(\mathrm{vN};\nu\geq
4)}(x)=\nu\left[-s_{\mathrm{vN}}(1)\frac{x}{2}+s_{\mathrm{vN}}(\tfrac{1}{\eta})\left(\frac{x}{4}\right)^{\tfrac{1}{\eta}}+s_{\mathrm{vN}}(4\eta)\left(\frac{x}{4}\right)^{4\eta}\right]+o\big{(}x^{\min(1,\tfrac{1}{\eta},4\eta)}\big{)}$
(5.16)

for the von Neumann tripartite information, where the function
$s_{\mathrm{vN}}$ is defined in Eq. (2.8). We note that, any time the
tripartite information is negative, the leading term is linear in $x$, as in
the non-interacting case $\eta=1/2$, also discussed in Ref. Fagotti2012New .

##  6 Spin structure from the fusion of models

#### With one-site shift invariance.

The tripartite information of the JW deformation has the form

$e^{(\alpha-1)[I_{3}^{(\alpha;\nu)}(x)+\log
2]}=\sum_{n=1}^{M_{\alpha}}g_{\alpha,n}\gamma_{\alpha,n}^{\nu}-\sum_{n=1}^{\tilde{M}_{\alpha}}\tilde{g}_{\alpha,n}(-\gamma_{\alpha,n})^{\nu}\,,$
(6.1)

where each $\gamma$ is a ratio of partition functions and $g_{\alpha,n}$ are
positive integers that represent the degeneracy of the terms. If we consider a
large enough number of deformations, we can generally express
$\\{\gamma_{\alpha,n}\\}$ in terms of
$\\{e^{(\alpha-1)[I_{3}^{(\alpha;\nu)}(x)+\log 2]}\\}$ by solving polynomial
equations. We have performed this procedure for $\alpha=2$ and $\alpha=3$ in
Appendix E. The result for $\alpha=2$ is given in Eq. (2.15).

#### With two-site shift invariance.

While a priori effective, the method with JW spin deformations has a rather
serious downside: it requires the analysis of models with next-to-next-to-
next-to-nearest-neighbor interactions. This can be numerically challenging
both because higher range is generally associated with higher numerical
complexity and because the typical lengths of a JW $\nu$-deformation is $\nu$
times the typical length of the original model. We can ease these problems by
generalising the deformations to the fusion of different models, as also
considered in Ref. Fagotti2012New . We can then set an upper bound to the
range, for example next-to-nearest-neighbor interactions, and consider a large
enough number of combinations of models that allow us to invert the problem.

Specifically, we consider the interleaved XXZ-XY model

$\bm{\mathbf{H}}_{\mathrm{XXZ}\oplus\mathrm{XY}}=\sum_{\ell}(1+\gamma)\bm{\mathbf{\sigma}}_{2\ell}^{x}\bm{\mathbf{\sigma}}_{2\ell+1}^{z}\bm{\mathbf{\sigma}}_{2\ell+2}^{x}+(1-\gamma)\bm{\mathbf{\sigma}}_{2\ell}^{y}\bm{\mathbf{\sigma}}_{2\ell+1}^{z}\bm{\mathbf{\sigma}}_{2\ell+2}^{y}+\bm{\mathbf{\sigma}}_{2\ell-1}^{x}\bm{\mathbf{\sigma}}_{2\ell}^{z}\bm{\mathbf{\sigma}}_{2\ell+1}^{x}+\bm{\mathbf{\sigma}}_{2\ell-1}^{y}\bm{\mathbf{\sigma}}_{2\ell}^{z}\bm{\mathbf{\sigma}}_{2\ell+1}^{y}\\\
+2h\bm{\mathbf{\sigma}}_{2\ell}^{z}+\Delta\bm{\mathbf{\sigma}}_{2\ell-1}^{z}\bm{\mathbf{\sigma}}_{2\ell+1}^{z}\,.$
(6.2)

The model is useful because in terms of Jordan-Wigner fermions
$\bm{\mathbf{c}}^{\dagger}_{\ell}=\prod_{j<\ell}\bm{\mathbf{\sigma}}_{j}^{z}\,\frac{\bm{\mathbf{\sigma}}_{\ell}^{x}+i\bm{\mathbf{\sigma}}_{\ell}^{y}}{2}$
it is a sum of the fermionic XXZ chain (2.2) on the odd sublattice and the
fermionic version of the quantum XY chain (3.34), i.e. the Kitaev chain, on
the even sublattice. For simplicity we work here with the XXZ chain in zero
magnetic field, and we denote by $h$ the magnetic field corresponding to the
XY chain, but the same considerations would apply to the XXZ with the magnetic
field.

Since the model is a sum of fermionic Hamiltonians we have the factorization
of the function $\mathcal{F}_{f}$, defined in Eq. (3.19), of the interleaved
model into the ones of the two constituting models. Let us denote by
$X|_{j}=X\cap(2\mathbb{Z}+j)$ for $j=0,1$ the restriction of a set $X$ to even
and odd sublattice respectively. We have then

$\mathcal{F}_{f}^{\mathrm{XXZ}\oplus\mathrm{XY}}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](A,B,C)=\mathcal{F}_{f}^{\mathrm{XXZ}}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](A|_{1},B|_{1},C|_{1})\mathcal{F}_{f}^{\mathrm{XY}}\left[\begin{smallmatrix}\vec{\varepsilon}\\\
\vec{\delta}\end{smallmatrix}\right](A|_{0},B|_{0},C|_{0})\,,$ (6.3)

where the superscripts indicate the model in which $\mathcal{F}_{f}$ is
computed, namely, the XXZ spin chain (2.1) and the XY spin chain (3.34).

The relation (6.3) can be exploited to compute $\mathcal{F}_{f}$ in the
original model, i.e. the XXZ spin chain, for any configuration $A,B,C$. For
simplicity we focus on the Rényi index $\alpha=2$, corresponding to the torus.
To compute
$\mathcal{F}_{f}^{\mathrm{XXZ}}\left[\begin{smallmatrix}\varepsilon\\\
\delta\end{smallmatrix}\right](A,B,C)$ we can consider the tripartite
information in the XXZ chain and in the interleaved model corresponding to
three different values of $\gamma,h$ close to or at the circle
$\gamma^{2}+h^{2}=1$. One could be $\gamma=1,h=0$, for which
$\mathcal{F}_{f}^{\mathrm{XY}}\left[\begin{smallmatrix}\varepsilon\\\
\delta\end{smallmatrix}\right](A,B,C)=(-1)^{4\varepsilon\delta}$
irrespectively of the configuration and the others should be taken not too far
from the multicritical point $\gamma=0,h=1$ (and not too close either). This
will give us a system of linear equations whose solution gives the quantity of
interest, i.e.
$\mathcal{F}_{f}^{\mathrm{XXZ}}\left[\begin{smallmatrix}\varepsilon\\\
\delta\end{smallmatrix}\right](A,B,C)$ for all torus spin structures
$\varepsilon,\delta$. We want to use the parameters $h,\gamma$ close to the
circle $h^{2}+\gamma^{2}=1$ because only there the dependence of
$\mathcal{F}_{f}$ on the configuration is sufficiently complex to give us
distinct enough coefficients of the linear problem (see Sec. 3.3).

Let us use the shorthand notation $\mathcal{C}=(A,B,C)$ for the configuration
of interest. Let us denote by $\tilde{A},\tilde{B},\tilde{C}$ three adjacent
subsystems that are twice the size of $A,B,C$, i.e.
$\tilde{A}=\\{1,2,\ldots,2|A|\\}$,
$\tilde{B}=\\{2|A|+1,2|A|+2,\ldots,2|A|+2|B|\\}$,
$\tilde{C}=\\{2(|A|+|B|)+1,2(|A|+|B|)+2,\ldots,2(|A|+|B|+|C|)\\}$, and let us
denote the whole configuration by
$\tilde{\mathcal{C}}=(\tilde{A},\tilde{B},\tilde{C})$. The interleaved model
satisfies
$\mathcal{F}_{f}^{\mathrm{XXZ}\oplus\mathrm{XY}}\left[\begin{smallmatrix}\varepsilon\\\
\delta\end{smallmatrix}\right](\tilde{\mathcal{C}})=\mathcal{F}_{f}^{\mathrm{XY}}\left[\begin{smallmatrix}\varepsilon\\\
\delta\end{smallmatrix}\right](\mathcal{C})\mathcal{F}_{f}^{\mathrm{XXZ}}\left[\begin{smallmatrix}\varepsilon\\\
\delta\end{smallmatrix}\right](\mathcal{C})$. Considering the configuration
$\mathcal{C}$ in the XXZ spin chain and the configuration
$\tilde{\mathcal{C}}$ in interleaved model for $(\gamma_{1},h_{1})=(1,0)$ and
other two values $(\gamma_{j},h_{j})$, $j=1,2$, close to the circle
$\gamma^{2}+h^{2}=1$, we obtain the system of equations

$\frac{1}{2}\left(\begin{smallmatrix}1&1&1&-1\\\ 1&1&1&1\\\
\mathcal{F}_{f}^{\mathrm{XY}_{2}}\left[\begin{smallmatrix}0\\\
0\end{smallmatrix}\right](\mathcal{C})&\mathcal{F}_{f}^{\mathrm{XY}_{2}}\left[\begin{smallmatrix}0\\\
1/2\end{smallmatrix}\right](\mathcal{C})&\mathcal{F}_{f}^{\mathrm{XY}_{2}}\left[\begin{smallmatrix}1/2\\\
0\end{smallmatrix}\right](\mathcal{C})&-\mathcal{F}_{f}^{\mathrm{XY}_{2}}\left[\begin{smallmatrix}1/2\\\
1/2\end{smallmatrix}\right](\mathcal{C})\\\
\mathcal{F}_{f}^{\mathrm{XY}_{3}}\left[\begin{smallmatrix}0\\\
0\end{smallmatrix}\right](\mathcal{C})&\mathcal{F}_{f}^{\mathrm{XY}_{3}}\left[\begin{smallmatrix}0\\\
1/2\end{smallmatrix}\right](\mathcal{C})&\mathcal{F}_{f}^{\mathrm{XY}_{3}}\left[\begin{smallmatrix}1/2\\\
0\end{smallmatrix}\right](\mathcal{C})&-\mathcal{F}_{f}^{\mathrm{XY}_{3}}\left[\begin{smallmatrix}1/2\\\
1/2\end{smallmatrix}\right](\mathcal{C})\end{smallmatrix}\right)\left(\begin{smallmatrix}\mathcal{F}_{f}^{\mathrm{XXZ}}\left[\begin{smallmatrix}0\\\
0\end{smallmatrix}\right](\mathcal{C})\\\
\mathcal{F}_{f}^{\mathrm{XXZ}}\left[\begin{smallmatrix}0\\\
1/2\end{smallmatrix}\right](\mathcal{C})\\\
\mathcal{F}_{f}^{\mathrm{XXZ}}\left[\begin{smallmatrix}1/2\\\
0\end{smallmatrix}\right](\mathcal{C})\\\
\mathcal{F}_{f}^{\mathrm{XXZ}}\left[\begin{smallmatrix}1/2\\\
1/2\end{smallmatrix}\right](\mathcal{C})\end{smallmatrix}\right)=\left(\begin{smallmatrix}\exp\left[I_{3}^{\mathrm{XXZ}}(\mathcal{C})\right]\\\
\exp\left[I_{3}^{\mathrm{XXZ}\oplus\mathrm{XY}_{1}}(\tilde{\mathcal{C}})\right]\\\
\exp\left[I_{3}^{\mathrm{XXZ}\oplus\mathrm{XY}_{2}}(\tilde{\mathcal{C}})\right]\\\
\exp\left[I_{3}^{\mathrm{XXZ}\oplus\mathrm{XY}_{3}}(\tilde{\mathcal{C}})\right]\end{smallmatrix}\right)\,,$
(6.4)

where by $\mathrm{XY}_{j}$ we denote the quantum XY chain (3.34) with
$(\gamma,h)=(\gamma_{j},h_{j})$ and we have also omitted writing the Rényi
index $\alpha$ in the tripartite information $I_{3}^{(\alpha)}$. The solution
to this linear system of equations unfolds the contribution from each the spin
structure in XXZ.

In practice, we compute the tripartite information numerically using tensor
network algorithms. In the interleaved model we take the values of
$\gamma_{j},h_{j}$ for $j=2,3$ exactly on the circle $\gamma^{2}+h^{2}=1$. The
accuracy of the numerical simulations is however crucial since a small
discrepancy in the parameters results in huge discrepancies in the
$\mathcal{F}$s. Thus, instead of computing $\mathcal{F}_{f}^{\mathrm{XY}}$
directly using analytical formulas (3.37), we ease such sensitivity by
renormalising the parameters as follows. We assume that the state is close to
the ground state of the XY model for some effective $h$ and $\gamma$ and
compare the predictions for the expectation values of two local observables
with their numerical values. For instance, we use the identities

$\begin{split}&\braket{\bm{\mathbf{\sigma}}^{z}_{2\ell}}=-\frac{1}{\pi}\int\limits_{0}^{\pi}\frac{h-\cos
k}{\sqrt{(\gamma\sin k)^{2}+(h-\cos k)^{2}}}dk\,,\\\
&\braket{\bm{\mathbf{\sigma}}^{x}_{2\ell}\bm{\mathbf{\sigma}}^{z}_{2\ell+1}\bm{\mathbf{\sigma}}^{x}_{2\ell+2}-\bm{\mathbf{\sigma}}^{y}_{2\ell}\bm{\mathbf{\sigma}}^{z}_{2\ell+1}\bm{\mathbf{\sigma}}^{y}_{2\ell+2}}=-\frac{2}{\pi}\int\limits_{0}^{\pi}\frac{\gamma\sin^{2}(k)}{\sqrt{(\gamma\sin
k)^{2}+(h-\cos k)^{2}}}dk\,,\end{split}$ (6.5)

which are easily obtained considering that the interleaved model is equivalent
to an XY chain when restricted to even sites. We extract the left hand sides
of (6.5) from the numerical simulations and invert the system of integral
equations for $h$ and $\gamma$. Finally, we use such effective $\gamma,h$ for
computing $\mathcal{F}_{f}^{\mathrm{XY}}$. Appendix D demonstrates explicitly
the importance of this step by comparing some results obtained with and
without it. The procedure for computing $\mathcal{F}_{f}^{\mathrm{XY}}$ has
been developed in Ref. Fagotti2010disjoint  and discussed also in Refs
Maric2023Universality ; Maric2023universality2 . For completeness, it is also
reported in Appendix C.

###  6.1 Tensor network implementation

In this section we describe the tensor network algorithms that we used to
obtain numerical data for the tripartite information of three adjacent
subsystems in the spin chains of interest. Let $\ket{\psi}$ be the ground
state of a given Hamiltonian $H$. After fixing a local basis of the Hilbert
space associated to each spin, we can expand the state as follows

$\ket{\psi}=\sum_{i_{1},...,i_{L}=1}^{d}T_{i_{1},i_{2},...,i_{L}}\ket{i_{1}}\ket{i_{2}}...\ket{i_{L}}\,,$
(6.6)

where $L$ is the chain size, $d$ is the dimension of the local Hilbert space
($d=2$ in our case of spin-$\frac{1}{2}$ chains). The ground state minimizes
the expectation value of the energy, which corresponds to minimising
$\braket{\psi}{\bm{\mathbf{H}}}{\psi}$ over the $d^{L}$-dimensional space of
the coefficients $T_{i_{1},i_{2},...,i_{L}}$. This is a hard numerical problem
to solve for growing $L$, hence one generally resorts to variational
principles. Instead of searching for the minimum over all the states of the
Hilbert space, one restricts the minimisation to the states that have a
particular form, referred to as the variational ansatz. The minimum over this
subset of the Hilbert space is obviously an approximation from above to the
true minimum. How well the variational minimum approximates the actual minimum
depends on the quality of the variational ansatz. The most flexible and
successful variational ansatz for 1D problems are arguably matrix product
states (MPS), which are states that can be parametrized as

$\ket{\psi}=\sum_{i_{1},...,i_{L}=1}^{d}\sum_{\ell_{1},...,\ell_{L-1}=1}^{\chi}(M^{(1)})^{i_{1}}_{\ell_{1}}(M^{(2)})^{i_{2}}_{\ell_{1},\ell_{2}}(M^{(3)})^{i_{3}}_{\ell_{2},\ell_{3}}...(M_{L})^{i_{1}}_{\ell_{L-1}}\ket{i_{1}}\ket{i_{2}}...\ket{i_{L}}\,,$
(6.7)

where $\chi$ is known as bond dimension of the MPS, and $M^{(j)}$ is a tensor
with two indices for $i\in\\{1,L\\}$ and three indices otherwise (see e.g.
Orus_TNreview ). For a given bond dimension, our ansatz has $\sim Ld\chi^{2}$
variational parameters. Therefore the bond dimension controls the
“expressivity” of the ansatz: a larger bond dimension corresponds to a larger
subset of the Hilbert space in the minimization problem (note that, in
principle, a bond dimension that scales exponentially in system size allows
for an exact representation of any state).

We are interested in working in the thermodynamic limit under the assumption
of invariance under shifts of two lattice sites (this is a symmetry common to
all the interleaved models that we considered). Shift-invariance allows us to
encode the full variational ansatz in just two distinct tensors with three
indices each, which should be repeated infinitely many times in order to
reconstruct the full state. In particular, following iTEBD_first ;
iTEBD_second , our MPS variational ansatz will look like

$\ket{\psi}=\sum_{...,i_{-1},i_{0},i_{1},i_{2},...=1}^{d}\sum_{...,\ell_{-1},\ell_{0},\ell_{1},\ell_{2},...=1}^{\chi}...M^{i_{-1}}_{\ell_{-2},\ell_{-1}}N^{i_{0}}_{\ell_{-1},\ell_{0}}M^{i_{1}}_{\ell_{0},\ell_{1}}N^{i_{2}}_{\ell_{1},\ell_{2}}...\ket{i_{-1}}\ket{i_{0}}\ket{i_{1}}\ket{i_{2}}...\,,$
(6.8)

for a total of $2d\chi^{2}$ variational parameters, contained in the two
tensors $M^{i}_{\ell,\ell^{\prime}}$ and $N^{i}_{\ell,\ell^{\prime}}$. To
simplify the notations, drawings are usually employed to represent the various
quantities occurring in the computation. Tensors are represented by rectangles
having as many legs as the number of their indices, while a shared leg between
two tensors stands for a sum over the corresponding index. For example, a
drawing representing the building blocks $M$ and $N$ is reported in Fig. 8a,
and the state (6.8) can be represented as in Fig. 8b.

Note that there is some gauge freedom in the MPS representation of a state:
given an invertible matrix $X$, the two tensors $XM^{i}X^{-1}$ and
$XN^{i}X^{-1}$ describe the same state as $M^{i}$ and $N^{i}$. One can use
this freedom in the parametrization to recast the state in a form that is
particularly convenient for computations, known as canonical form of the MPS.
When working in the thermodynamic limit, the canonical form can be defined
such that the contraction of the semi-infinite tails of the MPS simplifies to
the identity, as shown diagrammatically in Fig. 8c (see e.g. iTEBD_second  for
more details about canonical forms). This allows one to express many important
objects in a very compact form; as an example, we report in Fig. 8d the tensor
network representation of a reduced density matrix.

Figure 8: Conventional drawings representing tensor networks. Tensors are
represented by rectangles having as many legs as the number of their indices.
Since we assume 2-site-shift invariance, the MPS is parametrized by two
distinct 3-index tensors. a) The building block of the MPS ansatz, which is
denoted by $M$ and $N$ in the main text. b) MPS for the full state, given by
the infinite repetition of the building block, and its adjoint, where the
asterisk stands for complex conjugation. c) Canonical form of the MPS,
allowing one to simplify the contraction of the infinite tails of the MPS. d)
Example of a reduced density matrix for a block $A$ containing $6$ spins; for
the sake of brevity, we omitted the sum over the basis, analogous to the sums
in b). Figure 9: Tensor networks used to compute the Rényi-2 entropy. a)
Definition of the tensor representing the transfer matrix (the definition can
be easily adapted to the case in which the first site is associated to the
tensor $N$ instead of the tensor $M$ and/or to the case of an odd number of
sites). b) Definition of the twist field tensor; it is used, together with the
transfer matrix, to get more compact representations of the tensor networks
that follow. c) Entropy for a subsystem $A$ containing 6 spins. d) Entropy for
a subsystem made by two disjoint blocks $A$ and $C$ at distance $|B|$.

Once the variational ansatz (6.8) is chosen, the next step is to develop an
algorithm to find the parameters that minimize the energy. We use the
variational uniform matrix product states (VUMPS) algorithm, formulated in
Ref. VUMPSfirst  to perform variational energy minimization in infinite 1D
quantum systems. Ref. VUMPSfirst  shows that the VUMPS algorithm generally
outperforms its more-famous predecessors, the infinite Denisy Matrix
Renormalization Group (iDMRG) iDMRG  and the infinite Time Evolving Block
Decimation (iTEBD) iTEBD_first , both in convergence speed and accuracy of the
approximation. It also seems to be a good choice for studying the critical
systems we are interested in. The algorithm comes with two main schemes of
convergence to the variational minimum: sequential and parallel. The original
paper VUMPSfirst  showed that these two approaches are typically equivalent in
terms of performance and our simulations are based on the sequential
algorithm. In practice, we use the Julia-language implementation of the
algorithm available in the work-in-progress package ITensorInfiniteMPS.jl,
based on the ITensors.jl library ITensor . The results that we show in this
work have been obtained setting the bond dimension to $\chi=140$ or
$\chi=180$.

Note that the true GS in the (critical) systems that we investigated cannot be
represented as an MPS: on the one hand the entanglement of a bipartition in
two extensive subsystems grows as the logarithm of the subsystems size, which
is divergent in the thermodynamic limit; on the other hand, for any state
represented as an MPS with finite bond dimension, the entanglement follows an
area law and is finite (see e.g. Orus_TNreview ). Thus, a finite bond
dimension introduced an effective length effective_length_MPS1 ;
effective_length_MPS2 ; Pollmann2009Theory , and only the observables with a
range sufficiently smaller than that can be described by the approximate
ground state state well enough. We have checked the accuracy of the
approximation by verifying that, at the scales that we consider, the numerical
data are not significantly affected by an increase in the bond dimension.

The final step is to compute the Rényi entropy by contracting the tensor
networks shown in Fig. 9. The computational cost of the contraction scales as
$\chi^{5}$, which is much better than the scaling of the same computation in
an MPS with periodic boundary conditions, shown to be $\chi^{9}$ in Ref.
Coser2014OnRenyi . The reason is that here we can simplify the contraction by
employing the canonical form of the MPS, which is not possible with periodic
boundary conditions. This is another advantage of working directly in the
thermodynamic limit. Note that the green and pink tensors used in Fig. 9 to
give a compact representation of the tensor network of interest can be
interpreted as powers of the transfer matrix and twist fields respectively, as
was done e.g. in Ref. Coser2014OnRenyi . We do not discuss them here and we
simply use them as a convenient way to represent our tensor network.

###  6.2 On the numerical data

Figure 10: Contributions of different spin structures $(\varepsilon,\delta)$
to the Rényi-$2$ tripartite information of three finite adjacent subsystems
$A,B,C$ in the infinite XXZ spin chain with $\Delta=-0.3$. Function
$\mathcal{F}_{f}\left[\begin{smallmatrix}\varepsilon\\\
\delta\end{smallmatrix}\right](A,B,C)$ (see Eq. (3.18) and Eq. (3.19)) is
shown versus the cross ratio $x=|A||C|/[(|A|+|B|)(|B|+|C|)]$ for all four
torus spin structures. The lengths $|A|=|C|$ are varied, while the length of
$|B|$ is set to $6$ (filled markers) or $10$ (open markers). The points are
obtained using the method of fusion of models, incorporating the data from the
tensor network simulations. The solid lines are the predictions (2.10),
whereas dashed lines are the corresponding contributions in the free case
($\Delta=0$).  Figure 11: The Rényi-$2$ tripartite information in the infinite
fermionic XXZ chain (2) with $\Delta=-0.3$ for finite adjacent subsystems
$A,B,C$. The lengths $|A|=|C|$ are varied, while the length of $|B|$ is set to
$6$ (filled markers) or $10$ (open markers). The points are obtained using the
method of fusion of models, incorporating the data from the tensor network
simulations. The solid line is the prediction (2.4), whereas the dashed line
is the Rényi-$2$ tripartite information in the fermionic XX chain
($\Delta=0$), which is zero for any $x$.

We have used the described procedure with the interleaved model (6.2) to
compute the functions $\mathcal{F}_{f}\left[\begin{smallmatrix}\varepsilon\\\
\delta\end{smallmatrix}\right](A,B,C)$ in the ground state of the XXZ spin
chain (2.1) for all spin structures $\varepsilon,\delta$ of the torus, for
different configurations $(A,B,C)$ and various values of $\Delta$. We discuss
here the data for $(\Delta,h)=(-0.3,0)$, for which we reached a satisfactory
numerical convergence and are exemplar of the general behaviour observed.

Figure 4 showed the data at fixed cross ratio (1.4) equal to either $x=1/4$ or
$x=1/2$. In Fig. 10, instead, the length of the central subsystem $B$ is fixed
to two different values and the other lengths are set to the same value,
varying over a range that is numerically accessible with high enough accuracy.
Finally, Fig. (11) shows a comparison between prediction and numerical data
for the fermionic Rényi-$2$ tripartite information.

In all the cases, the data are in good visual agreement with the predictions
(2.10), but we point out that, in order to check the agreement quantitatively,
one would need to work out the asymptotics of the finite-size corrections,
which could depend on the particular spin structure and would not necessarily
match the leading corrections expected in the tripartite information. We have
not addressed this issue, therefore we leave a real quantitative check to
future investigations. Even if our analysis is more qualitative than
quantitative, it provides anyway the first numerical check of its kind in an
interacting model. It also provides an indirect check of the Rényi-$2$
tripartite information in the ground state of the Jordan-Wigner-$\nu$
deformations (2.3) for any $\nu$.

##  7 Insights into the tripartite information after global quenches

One of the commonest protocols of global quench starts with preparing the
system in the ground state of a local, translationally invariant Hamiltonian,
which is then suddenly perturbed, for example, by changing a global magnetic
field. In generic situations the system is expected to locally thermalise, in
the sense that the expectation values of local observables will relax to
thermal values. An exception to that behaviour, which has attracted much
interest, is when the system has infinitely many (pseudo)local conservation
laws, which generally means that the system is integrable. The stationary
state can then be characterised by a so-called generalised Gibbs ensemble
(GGE) Rigol2007Relaxation ; Doyon2017Thermalization ; Vidmar2016Generalized ;
Essler2016Quench . A GGE shares many properties in common with thermal states,
for example, subsystems are generally described by reduced density matrices
with an extensive entropy. It can however also exhibit exceptional features.
Specifically, static correlations in thermal states of systems with local
interactions decay exponentially with the distances
Bluhm2022exponentialdecayof ; in a generalised Gibbs ensemble of a system with
local interactions, that is not always the case. Ref. Maric2022Universality
considered some settings in which one should expect also algebraically
decaying correlations at late times after a global quench in an integrable
model. It was observed observed that such exceptional property leaves clear
marks in the entanglement entropies of subsystems. In particular, Refs
Maric2022Universality ; Maric2023Universality ; Maric2023universality2
studied the tripartite information after global quenches in noninteracting
spin chains and identified contributions that are generally killed by the
quench and others that are instead only mildly affected by nonequilibrium time
evolution. As a result, the tripartite information of large adjacent
subsystems is not trivialised by the limit of infinite time and, remarkably,
remains a function of the cross ratio such as in conformal critical ground
states. It was also pointed out that the final result has a universal
property: when the tripartite information survives the quench, the limit
$x\rightarrow 1^{-}$ turns out to be nonzero and, in the large class of
noninteracting models investigated in Ref. Maric2023Universality , it
approaches a universal $-\log 2$.

The analogous problem in the presence of interactions deserves a separate
investigation. Here we simply compute what would be the late time tripartite
information in 1D lattice models if the observations of Ref.
Maric2022Universality  about the quench having the main effect of killing the
expectation values of operators including fermionic strings between the blocks
held true in more general models. Namely, we do not expect this property to
depend on the boson radius in the underlying field theory. This has the effect
of killing all terms (3.19) in the tripartite information with $s_{k}\neq 1$
for some $k$. These are the terms with $\vec{\varepsilon}\neq 0$. We will also
assume the quench to be small, so that we can presumably ignore the
modifications in the remaining contributions.

If we keep only the terms that we expect to be nonzero we find

$I_{3}^{(\alpha)}(x)=\frac{1}{\alpha-1}\log\left(\frac{1}{2^{\alpha-1}}\sum\limits_{\vec{\delta}}\mathcal{F}_{f}\left[\begin{smallmatrix}\vec{0}\\\
\vec{\delta}\end{smallmatrix}\right](x)\right)=\frac{1}{\alpha-1}\log\frac{\Theta(\vec{0}|\tfrac{1}{\eta}\Omega)\Theta(\vec{0}|4\eta\Omega)}{[\Theta(\vec{0}|\Omega)]^{2}}\,,$
(7.1)

which can also be written as

$I_{3}^{(\alpha)}(x)=I_{3,\mathrm{GS}}^{(\alpha)}(x)+\frac{1}{\alpha-1}\log\frac{\Theta(\vec{0}|4\eta\Omega)}{\Theta(\vec{0}|\eta\Omega)}\,,$
(7.2)

where $I_{3,\mathrm{GS}}^{(\alpha)}$ is the tripartite information in the
ground state of the XXZ chain, given by Eq. (4.5).

#### The large $x$ expansion: residual tripartite information.

The large $x$ expansion for the tripartite information in the stationary state
can be obtained from the corresponding expansion of two terms in Eq. (7.2).
The large $x$ expansion of the first term, describing the ground state,
follows directly from the crossing symmetry
$I_{3,\mathrm{GS}}^{(\alpha)}(x)=I_{3,\mathrm{GS}}^{(\alpha)}(1-x)$ and the
small $x$ expansion, given by Eq. (4.6). It remains to deal with the large $x$
expansion of the second term in (7.2).

We start from the large $x$ expansion of $\beta_{k/\alpha}(x)$ appearing in
the expression (2.5) for the period matrix. The expansion follows from the
relation $\beta_{k/\alpha}(x)=1/\beta_{k/\alpha}(1-x)$ and the small $x$
expansion of $\beta_{k/\alpha}(x)$, reported in Refs
Calabrese2009Entanglement1 ; Calabrese2011Entanglement  and in Eq. (F.2) in
the appendix. Namely, we have for $x\to 1$

$\begin{split}&\beta_{\frac{k}{\alpha}}(x)=-\frac{\pi}{\sin\left(\pi\frac{k}{\alpha}\right)}\left[\log(1-x)+f_{\frac{k}{\alpha}}+\sum_{\ell=1}^{\infty}p_{\ell}\left(\frac{k}{\alpha}\right)(1-x)^{\ell})\right]^{-1},\\\
&f_{\frac{k}{\alpha}}\equiv
2\gamma_{E}+\psi\left(\frac{k}{\alpha}\right)+\psi\left(1-\frac{k}{\alpha}\right),\end{split}$
(7.3)

where $\gamma_{E}$ is the Euler gamma,
$\psi(z)\equiv\Gamma^{\prime}(z)/\Gamma(z)$ is the polygamma function and
$p_{l}(z)$ is a polynomial of degree 2$l$, whose explicit expression is not
needed. Factoring out the leading term $\log(1-x)$ and exploiting the formula
for the geometric series we obtain further

$\beta_{\frac{k}{\alpha}}(x)=\frac{\pi}{\sin\left(\pi\frac{k}{\alpha}\right)}\frac{1}{|\log(1-x)|}\left[1-\frac{f_{\frac{k}{\alpha}}}{|\log(1-x)|}+O\left(\frac{1-x}{|\log(1-x)|}\right)\right],\quad\textrm{as
}x\to 1\;.$ (7.4)

This gives the large $x$ expansion of the period matrix (2.5)

$\Omega=\frac{2\pi
i}{\alpha}\frac{\alpha\mathbb{I}-\vec{u}\vec{u}^{\mathrm{T}}}{|\log(1-x)|}\left[\mathbb{I}-\frac{\mathbb{I}+\vec{u}\vec{u}^{\mathrm{T}}}{\alpha|\log(1-x)|}B+O\left(\frac{1-x}{|\log(1-x)|}\right)\right],$
(7.5)

where $\vec{u}=(1,1,\ldots,1)^{\mathrm{T}}$ is a vector with $\alpha-1$