Entanglement content of kink excitations
Luca Capizzi Michele Mazzoni


 








 





equationsection



bibliography.bib

  


    ł()̊⟨⟩#1  #1
    #1  #1
    Entanglement content of kink excitations
    
        Luca Capizzi^1, Michele Mazzoni^2
        
    ==========================================
    
    ^1Université Paris-Saclay, CNRS, LPTMS, 91405, Orsay, France.
    ^2Department of Mathematics, City, University of London, 10 Northampton Square EC1V 0HB, UK.
    § ABSTRACT 
    
    Quantum one-dimensional systems in their ordered phase admit kinks as elementary excitations above their symmetry-broken vacua. While the scattering properties of the kinks resemble those of quasiparticles, they have distinct locality features that are manifest in their entanglement content. In this work, we study the entanglement entropy of kink excitations. We first present detailed calculations for specific states of a spin-1/2 chain to highlight the salient features of these excitations. Second, we provide a field-theoretic framework based on the algebraic relations between the twist fields and the semilocal fields associated with the excitations, and we compute the Rényi entropies in this framework. We obtain universal predictions for the entropy difference between the excited states with a finite number of kinks and the symmetry-broken ground states, which do not depend on the microscopic details of the model in the limit of large regions. Finally, we discuss some consequences of the Kramers-Wannier duality, which relates the ordered and disordered phases of the Ising model, and we explain why, counterintuitively, no explicit relations between those phases are found at the level of entanglement.
     
    
    § INTRODUCTION
    
    
    In the last decades, several studies regarding the entanglement of quantum systems shed light on the fundamental properties of quantum correlations in the ground state of many-body systems <cit.>. In particular, many exact results have been provided in the context of critical ground states via Conformal Field Theory (CFT) techniques <cit.>. Moreover, analytical investigations of gapped systems close to criticality have been given in a series of works by means of form-factor bootstrap for integrable systems <cit.>. Despite the abundance of works for quantum systems with a unique vacuum, for instance, many-body models in their disordered phase, only a few results are present for the symmetry-broken ground states, especially in the context of Quantum Field Theory (QFT). Considering the importance of spontaneous symmetry breaking in modern physics, this fact is quite surprising.
    
     Since there are Kramers-Wannier dualities connecting ordered and disordered phases <cit.>, one could erroneously argue that there should be no discernible distinction in the entanglement content between these phases, and consequently identical predictions should arise at corresponding dual points. However, this is not the case, as previous exact lattice results for the Ising chain pointed out <cit.>. While the presence of these discrepancies is known (see also Ref. <cit.>), some fundamental issues regarding their origin are not completely understood: is there a relation between the ground state entanglement of ordered and disordered phases? Is there a quantity dual to the entanglement entropy under Kramers-Wannier?
    
    In a previous work <cit.>, we have shown how the universal corrections to the area law found in Ref. <cit.> for a large class of theories with a single vacuum do not apply to the Ising Field Theory in its ordered phase (where spontaneous symmetry breaking arises and two vacua are present); that was traced back to the topological nature of the kink excitations and the corresponding selection rules for the form factors. A natural question, that is the focus of this work, is to understand whether similar discrepancies can be found for the low-lying excited states as well. Namely, in the series of works <cit.> the excess of entropy for the particle excitations in the disordered phase of a 1+1D gapped many-body system was proved to be universal. Those predictions boil down to general formulae with a simple semiclassical interpretation, and we ask whether similar results can be found for the ordered phase. We anticipate that we will find universal results different from those of the disordered phase; discrepancies are manifest, for example, in the Rényi entropies of tripartite geometries in the presence of a single kink.
    
    It is worth mentioning that the importance of the intrinsic non-locality of the kinks has been remarked for the quench dynamics in quantum spin chains in Refs. <cit.>. Moreover, topological frustration, as for antiferromagnetic chains with an odd number of sites, forces the ground state to be a one-kink state, with many counterintuitive consequences highlighted in Refs. <cit.>. Therefore, a deeper understanding of the entanglement content of kinks is imperative to correctly interpret a plethora of phenomena arising in one-dimensional quantum systems.
    
    We organize the manuscript as follows. In Section <ref> we first introduce the problem and obtain some explicit results on the lattice, and then, in the spirit of <cit.>, we provide a qubit picture that is sufficient to infer the correct universal contribution to the entanglement entropy arising from the kinks. Section <ref>, which is the core of our work, contains a field-theoretic explanation of the mechanism. In particular, we make use of the replica trick to compute the Rényi entropies via the so-called twist fields, and we trace back the origin of the novel entanglement properties to the algebraic relations between the twist fields and the disorder lines associated with the kinks. Finally, in Section <ref> we analyze the Kramers-Wannier duality and we explain the faith of the density matrix associated with local regions under the duality within a C^*-algebraic approach; this allows us to understand in detail why the entropy is not self-dual and we give an explicit construction of its dual. Conclusions are given in Section <ref>.
    
    
    
    
    § PREDICTIONS FROM SPIN CHAINS AND THE QUBIT PICTURE
    
    
    In this section, we compute the Rényi entropies of states of a 1D spin chain containing one or two kink excitations and we show that, in the case of a spatial tripartition, our predictions differ from those obtained for localised quasiparticle excitations. For the latter, the excess of Rényi and entanglement entropy with respect to the ground state is well known in the case of CFTs <cit.>, massive integrable QFTs <cit.> and spin chains <cit.>. In particular, in a state of a massive QFT containing a single localised excitation, the excess of entanglement entropy of a region A with fixed relative volume r= V_A / V in the large volume limit is given by 
    Δ S = -r log r - (1-r)log(1-r).
    
    This result relies on the locality of the excitation, and it also applies whenever A is made of disjoint connected regions. On the other hand, if we consider the topological excitations of a spin chain, such as the kinks in the ordered phase of the Ising model, the situation is different. Indeed, these excitations separate regions of different magnetisations, and thus by measuring the local magnetisation at a given point one could tell, in principle, whether the domain wall is found to the right or to the left of that point. The same is not possible for (local) quasiparticles, and physically one expects that this is the main feature responsible for the difference between kinks and particles at the level of entanglement content: in particular, we will show that (<ref>) does not hold for kinks if A consists of multiple disjoint regions.
    
    The situation becomes more complex in the presence of two or more kinks, leading to a richer phenomenology. In particular, if two kinks are sufficiently close they behave as a single collective particle (bound state): measuring the magnetisation at a distant point does not give information on the position of the bound state. Conversely, if the kinks are well separated, a measurement at a given point can reveal whether the two kinks belong to the opposite sides of the system. We depict the two scenarios in Figure <ref>.
    
    
    
    In this section, we study the low-lying spectrum of a spin chain deep in its order phase (e.g. the quantum Ising chain at small transverse field), where the two vacua become product states with opposite magnetisation and the kinks are plane-wave superpositions of domain wall configurations: this limiting regime, where the correlation functions of the vacua shrink to zero, is paradigmatic in its simplicity and it allows to elucidate the difference between quasiparticles and kinks. For states with one or two kinks, we are able to perform exact calculations on the lattice. However, since the task becomes challenging as the number of kinks increases, we also propose an alternative picture based on multi-qubit states; we adapt the qubit picture of Refs. <cit.>, initially introduced to describe the entanglement entropy and negativity of localised excitations, to deal with kinks.
    
    
    
     §.§ Rényi entropies of some one- and two-kink states in a spin-1/2 chain
    
    
      §.§.§ Structure of the Hilbert space
    
    
    We consider a spin-1/2 chain of length L, whose Hilbert space is
    ℋ = Span{|s_1,…,s_L⟩ ,    s_j = ±}≃łℂ^2^̊⊗ L.
    
    We regard the two states |+… +⟩ and |-… -⟩ as vacua, related by a global ℤ_2 symmetry, and we parametrise the other configurations via domain walls at different lattice sites as explained below. 
    
    First, according to the sign of the spin on the left boundary, we identify two sectors (of dimension 2^L-1) defined as
    ℋ^(±)≡Span{|±,s_2,…,s_L⟩,    s_j = ±},
    
    and such that ℋ = ℋ^(+)⊕ℋ^(-). Additionally, we split ℋ^(±) into sectors with a fixed number of (N) domain walls as
    ℋ^(±) = ⊕^L-1_N=0ℋ^(±)_N.
    
    Specifically, ℋ^(±)_N is generated by the configurations starting from ± with N changes of sign (domain walls) between neighboring sites.
    
    In the following, without loss of generality, we will focus on the sector ℋ^(+), with the first spin being up. For N=0, ℋ^(±)_0 is generated by |+,…,+⟩, meaning that no domain walls are present. For N=1, ℋ^(+)_1 admits the following basis
    |K_+-(j)⟩≡|+,+,…,+,-,…,-⟩,    j=1,…,L-1,
    
    with a domain wall between the positions j and j+1. Similarly, for N=2, ℋ^(+)_2 is spanned by the states
    |K_+-(j_1)K_-+(j_2)⟩≡|+,+…,+,-,…,-,+,…,+⟩,     1 ≤ j_1 <j_2 ≤ L-1,
    
    with the first domain wall between sites j_1, j_1 +1, and the second between sites j_2 and j_2 +1. More generally, a basis for ℋ^(+)_N is parametrised by sequences of N ordered positions (j_1,…,j_N), belonging to {1,…,L-1}, and therefore the dimension of this sector is
    dimℋ^(+)_N = L-1N;
    
    this is consistent with ∑^L-1_N=0dimℋ^(+)_N = dimℋ^(+), since ∑_N=0^L-1L-1N=2^L-1.
    
    
    
      §.§.§ One-kink state
     
    
    In this section, we analyse a state given by a coherent superposition of configurations with a single domain wall:
    |𝒦_+-(p)⟩ = 1/√(L-1)∑^L-1_j=1e^i p j|K_+-(j)⟩.
    
    This state represents a kink with a given momentum p that is completely delocalised in space. 
    
    We are interested in the computation of the Rényi entropy of extended regions, and we first analyse the simplest case of a bipartition A∪A̅, with
    
    A = {1,…,ℓ},  A̅ = {ℓ +1,… L}.
    
    We can compute the reduced density matrix (RDM) by tracing out the degrees of freedom in A̅:
    ρ^(1)_A = Tr_A̅ł|𝒦_+-(p)⟩⟨𝒦_+-(p)|≡̊s_ℓ + 1,…,s_L = ±∑A̅⟨s_ℓ +1,…,s_L | 𝒦_+-(p)|⟩⟨𝒦_+-(p)| s_ℓ +1,…,s_L|_⟩A̅,
    
    obtaining
    
        ρ^(1)_A    = 1/L-1∑_j,j'=1^ℓ-1 e^i p (j-j')|K_+-(j)⟩_A A⟨K_+-(j')| + ł1-ℓ/L-1|̊%̊s̊⟩̊+,…,+_A A⟨+,…,+|
           +1/L-1∑_j=1^ℓ-1ł e^i p (j-ℓ)|K_+-(j)⟩_A A⟨+,…,+| + h.c.,̊
    
    with |K_+-(j)⟩_A the state of A with a kink at position j. We identify three types of contributions in Eq. (<ref>): a term with one kink in A, a term with no kinks in A, and a mixed term that couples the two sectors; in particular, one can show that the latter contribution is due to the configurations with a domain wall lying between ℓ and ℓ+1, that is the edge of the region A. While the exact diagonalisation of ρ^(1)_A at finite size is quite cumbersome, a simplification occurs in the limit
    ℓ→∞,    L →∞,     r≡ℓ/L  fixed.
    
    In particular, the first two terms of (<ref>) can be diagonalised simultaneously and the associated non-vanishing eigenvalues are r and 1-r respectively: the third term can be treated in perturbation theory, and in the limit (<ref>) one can easily show that the associated contribution to the spectrum goes to zero. Therefore, in this limit, the nth Rényi entropy of the one-kink state is
    
    S_n = 1/1-nlog[(ρ_A^(1))^n] ≃1/1-nlog(r^n + (1-r)^n).
    
    We observe that this result coincides with the predictions of Ref. <cit.> for quasiparticles; in particular, the entanglement entropy, obtained in the limit n→ 1 of Eq. (<ref>), directly yields (<ref>) since the entropy of the vacuum state |+,…,+⟩ vanishes. This means that in the case of the geometry (<ref>) there is no distinction between kinks and quasiparticle at the level of entanglement entropy (at least, in the large volume limit).
    
    We now turn to a tripartite geometry, and we consider the entanglement between a single interval and the rest. Namely, we choose A = A_1∪ A_3, A̅=A_2 with
    
    A_1 = {1,…,ℓ_1},    A_2 = {ℓ_1 +1,…,ℓ_1 + ℓ_2},    A_3 = {ℓ_1 + ℓ_2 +1,…,L};
    
    therefore, after defining ℓ_3 ≡ L-ℓ_1-ℓ_2, the length of the interval A_i is ℓ_i.
    The evaluation of the reduced density matrix ρ^(1)_A_1∪ A_3 is analogous to the bipartite case discussed above. Thus, after some algebra, we obtain the following expression
    
        ρ^(1)_A_1 ∪ A_3 = 1/L-1∑_j,j'=1^ℓ_1-1 e^ip (j-j')|K_+-(j)⟩_A_1A_1⟨K_+-(j')|⊗|-…,-⟩_A_3A_3⟨-,…,-|
           + 1/L-1∑_j,j'=ℓ_1+ℓ_2 +1^L-1 e^ip(j-j')|+,…,+⟩_A_1A_1⟨+,…,+|⊗|K_+-(j)⟩_A_3A_3⟨K_+-(j')|
           + ℓ_2 +1/L-1|+,…,+⟩_A_1A_1⟨+,…,+|⊗|-…,-⟩_A_3A_3⟨-,…,-|
           + 1/L-1∑_j=1^ℓ_1 -1ł e^i p (j-ℓ_1)|K_+-(j)⟩_A_1A_1⟨+,…,+|⊗|-…,-⟩_A_3A_3⟨-,…,-| + h.c.
           + 1/L-1∑_j=ℓ_1 + ℓ_2 + 1^L-1ł e^i p (j-ℓ_1-ℓ_2)|+,…,+⟩_A_1A_1⟨+,…,+|⊗|-…,-⟩_A_3A_3⟨K_+-(j)| + h.c..̊
    
    As expected, the number of kinks belonging to A_1 ∪ A_3 can be either 0 or 1, and we identify the first three terms in (<ref>) as those associated with a fixed number of kinks in the regions A_1 and A_3; specifically, after decomposing the Hilbert space ℋ_A_1⊗ℋ_A_3 in sectors, we find
    
      *  One (ℓ_1 -1) × (ℓ_1-1) block acting on ℋ_A_1,1^(+)⊗ℋ_A_3,0^(-)  (one kink in A_1, no kinks in A_3),
        
      *  One (ℓ_3 -1) × (ℓ_3-1) block acting on ℋ_A_1,0^(+)⊗ℋ_A_3,1^(+)  (no kinks in A_1, one kink in A_3),
        
      *  One 1 × 1 block on acting ℋ_A_1,0^(+)⊗ℋ_A_3,0^(-) (no kinks in A).
    
    On the other hand, the last two terms of (<ref>) mix the sectors. However, in the limit L→∞ with
    
    r_i ≡ℓ_i/L fixed,    i=1,2,3;
    
    the eigenvalues of the first three blocks converge to r_1, r_2, r_3 respectively, and the contribution coming from the two mixing terms goes to zero. In conclusion, we compute the Rényi entropy as
    
    S_n = 1/1-nlog[(ρ_A_1∪ A_3^(1))^n] ≃1/1-nlog(r_1^n + r_2^n + r_3^n),
    
    which coincides with (<ref>) for r_3=0. As anticipated, this result differs from the one of quasiparticles in a tripartite geometry (<ref>), where no distinction between connected and disconnected regions is found. The difference between the two results is a hallmark of the non-locality of kink excitations, and it is physically better understood when one takes the limit r_2 → 0: we shall come back to this point at the end of Section <ref>.
    
    
    
      §.§.§ Two-kink states
    
    
     Interesting features related to the non-locality of kink excitations emerge in the presence of two kinks for the tripartite geometry in Eq. (<ref>). As a paradigmatic example, we focus on the case where the two kinks have equal momenta, p=p'=0, and we consider the state
     |𝒦_+-(0)𝒦_-+(0)⟩=√(2/(L-1)(L-2))∑_1≤ j < j' ≤ L-1|K_+-(j)K_-+(j')⟩,
     
    where the prefactor ensures the correct normalisation. The computation of its reduced density matrix, denoted by ρ^(2)_A_1 ∪ A_3, does not pose additional difficulties, but it is lengthy and requires some care. We leave the details of the computation and the explicit expression of ρ^(2)_A_1 ∪ A_3 to Appendix <ref>; here, we discuss its salient features.
    
    First, since two kinks are present in the entire system, only zero, one, or two kinks can be present in A_1 ∪ A_3. Second, in analogy to the previous calculation, only the contributions coming from the blocks with a fixed number of kinks in A matter in the limit (<ref>); we report them below, together with their non-zero eigenvalues:
    
      *  One 1 × 1 block acting on ℋ_A_1,0^(+)⊗ℋ_A_3,0^(+) (no kinks in A). It yields an eigenvalue r_2^2. 
        
      *  One (ℓ_1 -1) × (ℓ_1-1) block acting on ℋ_A_1,1^(+)⊗ℋ_A_3,0^(+)  (one kink in A_1, one in A̅). It yields an eigenvalue 2r_1 r_2.
        
      *  One (ℓ_3 -1) × (ℓ_3-1) block acting on ℋ_A_1,0^(+)⊗ℋ_A_3,1^(-)  (one kink in A̅, one in A_3). It yields an eigenvalue 2r_3 r_2.
        
      *  One (ℓ_1 -1)(ℓ_3-1) × (ℓ_1 -1)(ℓ_3-1) block acting on ℋ_A_1,1^(+)⊗ℋ_A_3,1^(-)  (one kink in A_1 and one in A_3). It yields an eigenvalue 2r_1 r_3.
        
      *  One [(ℓ_1 -1)(ℓ_1-2)/2 + (ℓ_3 -1)(ℓ_3-2)/2] ×[(ℓ_1 -1)(ℓ_1-2)/2 + (ℓ_3 -1)(ℓ_3-2)/2] block acting on[Here, the tensor product ⊗ is performed first, followed by the direct sum ⊕.]ℋ_A_1,2^(+)⊗ℋ_A_3,0^(+)⊕ℋ_A_1,0^(+)⊗ℋ_A_3,2^(+) (two kinks in A). It yields an eigenvalue r_1^2 + r_3^2.
    
    
    It is worth noting that the last block gives rise to a mixing of two sectors, corresponding to the two kinks being in either A_1 or A_3, respectively. From the point of view of entanglement, this contribution can be interpreted as the two kinks behaving like a single bound state which has a probability r_1^2 + r_3^2 of being found in A_1 ∪ A_3. Similarly, the first block corresponds to a single bound state belonging to A_2 with probability r^2_2. Conversely, in all the other cases the kinks belong to distinct regions, and their genuine non-locality manifests explicitly in their entanglement contributions.
    
    We conclude by reporting the value of the Rényi entropy for the tripartition in the large volume limit:
    
    S_n ≃1/1-nlog(r_2^2n + (r_1^2 + r_3^2)^n + 2^n(r_1^n r_2^n + r_1^n r_3^n+r_2^n r_3^n)).
    
    This expression is different from that of the two quasiparticle excitations in Ref. <cit.>, which depends only on the sum r_1 + r_3 and not on r_1 and r_3 separately.
    
    
    
    
     §.§ The qubit picture
    
    In this section we further explore the difference between multi-kink and multi-particle states at the level of entanglement by adopting the qubit formalism, first introduced in <cit.>, and later extended to deal with internal symmetries <cit.>. The idea is to recast the states considered in Section (<ref>), into simpler states by retaining the same entanglement content (at least, in the large volume limit). Before showing how to generalise the formalism to kinks, we first review the qubit description in the case of quasiparticle excitations.
    
    
    
      §.§.§ One- and two-particle states, revisited
    
    Let us consider the same tripartite geometry of Section (<ref>), and we focus on a single quasiparticle at a given momentum (completely delocalised in space). We assign 1 to the region A_i whenever the particle belongs to it, and 0 otherwise. The qubit representation of the state above is
    |Ψ^(1)⟩ = √(r_1)|100⟩ + √(r_2)|010⟩ + √(r_3)|001⟩,
    
    so that r_i is the probability that the particle belongs to A_i.
    Its RDM is
    ρ_A^(1) = _A_2|Ψ^(1)⟩⟨Ψ^(1)| =
    ł[                  |00⟩       |10⟩       |01⟩;       ⟨00|        r_2          0          0;       ⟨10|          0        r_1 √(r_1 r_3);       ⟨01|          0 √(r_1 r_3)        r_3;            ].̊
    After introducing
    |Φ_0⟩ = |00⟩,   |Φ_1⟩ = √(r_1)|10⟩+√(r_3)|01⟩/√(r_1 + r_3),
    
    we express the RDM in terms of its spectral projectors as
    ρ_A^(1) = (1-r)|Φ_0⟩⟨Φ_0| + r|Φ_1⟩⟨Φ_1|,
    
    with r=r_1+r_3. In particular, the entanglement content, encoded in the non-zero eigenvalues of ρ^(1)_A (r and 1-r) depends only on the probability r that the particle belongs to A=A_1∪ A_3.
    
    The situation is slightly more complex in the presence of two indistinguishable particles. In that case, we assign the value 2 to the region A_i if both particles belong to A_i and we express the qubit-state as
    |Ψ^(2)⟩ = r_1 |200⟩ + r_2|020⟩ + r_3 |002⟩ + √(2r_1 r_2)|110⟩ + √(2r_1 r_3)|101⟩+√(2r_2 r_3)|011⟩,
    
    whose corresponding RDM is
    
        ρ_A^(2)   = _A_2|Ψ^(2)⟩⟨Ψ^(2)|
           =
        ł[                              |00⟩             |10⟩             |01⟩             |20⟩             |02⟩             |11⟩;             ⟨00|            r_2^2                0                0                0                0                0;             ⟨10|                0        2 r_2 r_1  2 r_2√(r_1 r_3)                0                0                0;             ⟨01|                0 2 r_2 √(r_1 r_3)        2 r_2 r_3                0                0                0;             ⟨20|                0                0                0            r_1^2          r_1 r_3 r_1 √(2 r_1 r_3);             ⟨02|                0                0                0          r_1 r_3            r_3^2 r_3 √(2 r_1 r_3);             ⟨11|                0                0                0 r_1 √(2 r_1 r_3) r_3 √(2 r_1 r_3)        2 r_1 r_3;                  ].̊
    
    Its non-zero eigenvalues are (1-r)^2, 2r(1-r), r^2 and, again, they only depend on r_1,r_3 though their sum r.
    An important message is that mixing is always present between the states with a given number of particles in A_1∪ A_3: as we will see, this is not the case for the kinks.
    
    
    
      §.§.§ Rényi entropies of multi-kink states in a tripartite geometry
    
    
    We now adopt a qubit-like formalism to describe the states containing multiple kinks/antikinks. The key idea consists of fixing the number of kinks in a given region and then forgetting about all the other details related to the spatial distributions of spins in that region. To this aim, we introduce the states |±⟩ to represent the two vacua, and |1^±⟩ as a kink (or antikink) interpolating between |±⟩ and |∓⟩. The state of a single kink delocalised in space is
    |Ψ^(1)⟩ = √(r_1)|1^+, -, -⟩ +√(r_2)|+, 1^+, -⟩+ √(r_3)|+, +, 1^+⟩,
    
    and its RDM for the region A=A_1∪ A_3 is
    
        ρ_A^(1) = r_1 |1^+,-⟩⟨1^+,-| + r_2 |+,-⟩⟨+,-| + r_3 |+,1^+⟩⟨+,1^+|.
    
    Crucially, this result reproduces the same non-zero spectrum of (<ref>) in the large-volume limit and it shows the absence of mixing between the states with a single kink in A_1 or A_3 respectively (in contrast with the quasiparticle state in Eq. (<ref>)). 
    
    We now focus on a two-kink state. First, we introduce the state |2^±⟩ to represent the presence of a kink-antikink pair interpolating between |±⟩ and |±⟩; then we consider
    
        |Ψ^(2)⟩   = r_1 |2^+, +, +⟩ +r_2 |+, 2^+, +⟩+ r_3 |+, +, 2^+⟩
           + √(2 r_1 r_2)|1^+, 1^-, +⟩ + √(2 r_1 r_3)|1^+, -, 1^-⟩ + √(2 r_2 r_3)|+, 1^+, 1^-⟩,
    
    whose RDM is
    
        ρ_A^(2) =
        ł[               |+,+⟩   |+,1^-⟩   |1^+,+⟩   |2^+,+⟩   |+,2^+⟩ |1^+,1^-⟩;     ⟨+,+|     r_2^2         0         0         0         0         0;   ⟨+,1^-|         0 2 r_2 r_3         0         0         0         0;   ⟨1^+,+|         0         0 2 r_1 r_2         0         0         0;   ⟨2^+,+|         0         0         0     r_1^2   r_1 r_3         0;   ⟨+,2^+|         0         0         0   r_1 r_3     r_3^2         0; ⟨1^+,1^-|         0         0         0         0         0 2 r_1 r_3;           ].̊
    
    In analogy with Section <ref>, a mixing occurs between the state with a pair kink-antikink in A_1 and that with the pair in A_3. For completeness, we diagonalise the RDM as
    ρ_A^(2) = r_2^2 |Φ_0^+⟩⟨Φ_0^+| + 2 r_1 r_2 |Φ_1^-⟩⟨Φ_1^-| + 2 r_2 r_3 |Φ_1^+⟩⟨Φ_1^+| + 2 r_1 r_3 |Φ_2^-⟩⟨Φ_2^-| + (r_1^2 + r_3^2) |Φ_2^+⟩⟨Φ_2^+|,
    
    with |Φ_0^+⟩ = |+,+⟩, |Φ_1^-⟩ = |1^+,+⟩, |Φ_1^+⟩=|+,1^-⟩, |Φ_2^-⟩=|1^+,1^-⟩ and
    |Φ_2^+⟩ = r_1 |2^+,+⟩ + r_3 |+,2^+⟩/√(r_1^2 + r_3^2).
    
    We find, as expected, the same non-zero spectrum reported in Section <ref>.
    
    Finally, we discuss the situation of a generic number N of kinks. In analogy with the case of indistinguishable particles <cit.>, we consider an N-kink state
    
        |Ψ^(N)⟩   = ∑_N_1, N_2, N_3 ∈ℕ
         N_1 + N_2 + N_3 = N√(N! r_1^N_1r_2^N_2r_3^N_3/N_1! N_2! N_3!)δ_N_1 + N_2 + N_3, N|N_1^+, N_2^(-)^N_1, N_3^(-)^N_1+N_2⟩_Kinks in all the three regions
           + ∑_N_1=1^N-1√(NN_1r_1^N_1r_2^N-N_1)|N_1^+, (N-N_1)^(-)^N_1,(-)^N⟩
           + ∑_N_1=1^N-1√(NN_1r_1^N_1r_3^N-N_1)|N_1^+, (-)^N_1,(N-N_1)^(-)^N_1⟩
           + ∑_N_1=1^N-1√(NN_1r_2^N_1r_3^N-N_1)|+, N_1^+,(N-N_1)^(-)^N_1⟩_Kinks in two regions
           + r_1^N/2|N, (-)^N,(-)^N⟩+r_2^N/2|+, N^+, (-)^N⟩+r_3^N/2|+, +, N^+⟩_Kinks in one region.
    
    In the expression above |N^±_j⟩ is the state with N_j kinks in a given region with the vacuum |±⟩ on their left. The coefficients of the state are the square roots of the probabilities of finding a certain configuration of kinks among the regions A_1, A_2, A_3. While the RDM does not have a very insightful structure, we write its diagonal form as
    ρ_A^(N) = λ_0 |Φ_0⟩⟨Φ_0| + ∑_ϵ=±
     N_A=1,…,Nλ_N_A^ϵ|Φ_N_A^ϵ⟩⟨Φ_N_A^ϵ|.
    
    The eigenstates |Φ_N_A^ϵ⟩ have a fixed number N_A of kinks in A_1 ∪ A_3 and a fixed parity ϵ of the number of kinks in one of these two regions: without loss of generality, we choose A_1, so ϵ =+ iff there is an even number of kinks in A_1. The parity of the number of kinks/antikinks in A_3 is then automatically fixed by N_A and N. In the eigenspace with N_A and ϵ fixed, the qubit states are combined in a superposition to form |Φ_N_A^ϵ⟩, as in Eq. (<ref>). 
    The mechanism behind the mixing is the following: two states form a linear superposition if and only if they have the same number N_Aand they can be obtained one from the other by moving one magnon (i.e. a pair kink/antikink) from one region to the other. For every N_A  0 there are two eigenspaces labeled by the parity ϵ, while for N_A = 0 there is a single one-dimensional eigenspace. Specifically, the eigenvalues are:
    λ_0 = r_2^N,
    
    and for N_A=1,…,N
        λ_N_A^+    = ∑_j=0,
        j  even^ N_AN!/j!(N_A - j)! (N-N_A)!r_1^j r_3^N_A - jr_2^N-N_A,
        λ_N_A^-    = ∑_j=1, j  odd^ N_AN!/j!(N_A - j)! (N-N_A)!r_1^j r_3^N_A - jr_2^N-N_A,
    
    where in every sum j, N_A-j and N-N_A correspond to the number kinks in A_1, A_3, A_2 respectively. 
    From (<ref>) and (<ref>) we can finally express the tripartite nth Rényi entropy of the N-kink state as
    
    S_n = 1/1-nlog[λ_0^n + ∑_N_A = 1^N ł(λ_N_A^+)^n + (λ_N_A^-)^n].
    
    The above expressions can be checked via a direct calculation; in Appendix <ref>, as a detailed example, we discuss explicitly the cases with N=3,4.
    
    
    
    
    
    § FIELD-THEORETIC APPROACH
    
    
    
    In this section, we discuss the origin of the unusual entanglement content of kink excitations within a general field-theoretical framework. Our discussion relies on the description of those excitations in terms of semilocal fields acting on a symmetry-broken ground state, and their relations with the twist fields, which are the building blocks to compute entanglement measures via replica trick. In particular, we characterise the algebraic relations between the twist fields and semilocal fields such as the disorder operator in the Ising QFT. We recall that, as explained in Ref. <cit.>, the algebraic relations between twist fields and local operators are sufficient to recover the entanglement content of quasi-particles; here, by generalising the commutation relations so to take into account the semilocality of the disorder operator, we are able to obtain the entanglement content of kink excitations.
    
    We consider a 1+1 QFT displaying spontaneous symmetry breaking of a ℤ_2 symmetry, for instance the Ising model in the ferromagnetic phase, and we denote its vacua by |±⟩. We take a finite-size system of length L, and we specify the boundary conditions at the extremal points x=0,L. In particular, we denote by b_± the boundary conditions associated with the presence of positive/negative boundary magnetic field. We assume that the asymptotic spectrum of the theory is described by kinks interpolating between the two vacua. Due to topological constraints, if the boundary conditions at x=0,L are equal, say they are both b_+, the ground state is the vacuum |+⟩ and kinks above it can only appear in an even number. On the other hand, if the system has different boundary conditions at the two edges, it can only host an odd number of kinks and the ground state of the finite-size system is a one-kink state.
    
    Let μ(x) be a disorder field, that is a semilocal field connecting the two vacua in the region x'>x. Roughly speaking [
    In the context of Integrable Field Theories <cit.>, one usually refers to μ as a field with UV scaling dimension Δ_μ = 1/8, while the elementary excitations correspond to the fermionic field with dimension Δ = 1/2; there is no direct relation between those fields, and μ has non-vanishing form factors with an arbitrary odd number of kinks. However, only the semilocal properties of the fields will enter our calculation and these issues do not play a role.], μ(x) generates a kink at position x, and we are interested in the associated state at a given momentum. A single kink state can only connect sectors with different boundary conditions. For instance, the state
    μ(x)|+⟩,
    
     interpolates between b_+ (at x=0) and b_- (at x=L). We now construct the Fourier transform of μ(x) as
    μ(p) =  ∫ dx  (e^-ipx +ℛe^ipx)μ(x),
    
    which generates a kink excitation at a given momentum (in absolute value). Here, ℛ is a phase (|ℛ|=1) which depends on the boundary conditions of the field μ, but it does not play a major role in the discussion below. We describe eigenstates with a finite number of kinks (in the large volume limit) as
    μ(p_1)…μ(p_2N)|+⟩,     (b_+,b_+)  sector,
    μ(p_1)…μ(p_2N)|-⟩,     (b_-,b_-)  sector,
    μ(p_1)…μ(p_2N+1)|+⟩,     (b_+,b_-)  sector,
    μ(p_1)…μ(p_2N+1)|-⟩,     (b_-,b_+)  sector.
    
    We remark that, in principle, the momenta associated with eigenstates of the finite-size systems have to be quantised; however, this issue does not play a role as long as the number of particles is finite and the momenta are kept fixed in the large volume limit (as in Ref. <cit.>).
    
    Let us now define smeared fields by restricting the support of μ(p) to some spatial region A. That is, we introduce
    μ_A(p) =  ∫_A dx  (e^-ipx +ℛe^ipx)μ(x).
    
    Then, we study the product of two smeared fields in the semiclassical limit where the momenta are fixed and the size of the regions becomes very large:
    łμ_A'(p') ^̊†μ_A(p) = ∫_A dx ∫_A' dx'  (e^-ipx +ℛe^ipx)  (e^ip'x' +ℛe^-ip'x') μ^†(x') μ(x) ≃
    ∫_A dx ∫_A' dx'  (e^-ipx +ℛe^ipx)  (e^ip'x' +ℛe^-ip'x') μ^†(x')μ(x), 
    
    where we only kept the lightest contribution coming from the fusion μ^†×μ→ 1, as all terms corresponding to heavier operators are less relevant in the limit of large regions; thus, we replace μ^†(x') μ(x) with its vacuum expectation value over |+⟩ (a similar discussion can be found in Ref. <cit.> with periodic boundary conditions for local fields). We can further simplify Eq. (<ref>) if we assume that μ^†(x)μ(x') decays fast enough with the distance |x-x'|, e.g. exponentially as it happens in the massive (ordered) phase of the Ising model. In that case, we can perform the change of variable x' = x+x” and approximate
    łμ_A'(p') ^̊†μ_A(p) ≃∫ dx” ∫_A ∩ A'dx (e^i(p'-p)x+ip'x” + ℛe^i(p+p')x+ip'x”+c.c.)μ^†(x)μ(x+x”)≃
    ∫ dx”ł V_A ∩ A'δ_p,p'e^ipx” + V_A ∩ A'δ_p+p',0ℛ + c.c.μ^†(x)μ(x+x”),
    
    where the vacuum correlation function μ^†(x)μ(x+x”) is assumed to be x-independent if x is sufficiently distant from the boundary points and we denoted by V_A ∩ A' the volume of the region A ∩ A'. Eventually, we get
    łμ_A'(p') ^̊†μ_A(p) ∝ V_A ∩ A'δ_p,p',
    
    where the proportionality constant comes from the integration of the vacuum correlation function. 
    
    We remark that the considerations above hold only if the momentum is kept fixed in the large-volume limit. However, in some relevant situations, this is not the case. For example, the lowest-lying one-kink state is described by a single-body wave function ∼sinłπ x/L$̊ (p=π/Landℛ=-1in Eq. (<ref>)) if Dirichlet boundary conditions are present. In particular, the density of particles is not uniform in the large-volume limit, and the state looks inhomogeneous. These cases can be tackled as well with minor modifications. In particular, one can consider general smearing induced by a functionf(x)as
    μ_f ≡∫ dx  f(x)μ(x).
    
    Inhomogeneous states can then be studied by computing correlation functions of these smeared fields; however, this is beyond the purpose of our work.
    
    
    
    
     §.§ Twist fields and their algebra
    
    
    We first review some properties of the twist fields and their relations with the local operators, following Refs. <cit.>, and then we generalise the discussion to the case of the semilocal disorder operator. We anticipate that the twist fields relevant to our discussion are those with an even number of disorder lines attached to them and associated with distinct copies of the replica field theory. To the best of our knowledge, these fields have been studied for the first time in our previous work <cit.> in the context ofℤ_2entanglement asymmetry; in particular, they have the same scaling dimension of the standard twist field but different monodromy properties, which become manifest in the ordered phase. We also mention that in previous works, for the symmetry-resolved entanglement, twist fields attached to a single disorder line have been studied in the disordered phase <cit.>: these are the composite branch-point twist fields obtained from the fusion of the standard field with the disorder operator.
    
    Let us considernreplicas of the QFT. A field𝒪(x)inserted in thejth replica is denoted by𝒪^j(x). We focus on an algebra𝒜of local fields; in the Ising model, this is the algebra generated by the order operatorσ(x)and the energy densityε(x). To compute the entropy of a spatial region, that is the entropy of the subalgebra of𝒜associated with that region, we employ the branch-point twist fields. These are fields in the replica model satisfying
    𝒯(x)𝒪^j(y) = 𝒪^j(y)𝒯(x),    x>y,
    𝒯(x)𝒪^j(y) = 𝒪^j+1(y)𝒯(x),    x<y.
    
    Loosely speaking, we say that𝒯(x)introduces a branch cut atx'>xwhich connects the replicas through the permutationj→ j+1. Similarly, we introduce a conjugate twist field𝒯̃(x)such that
    𝒯̃(x)𝒪^j(y) = 𝒪^j(y)𝒯̃(x),    x>y,
    𝒯̃(x)𝒪^j(y) = 𝒪^j-1(y)𝒯̃(x),    x<y
    
    acting as the permutationj→ j-1onx'>x.𝒯and𝒯̃are the building blocks to construct the Rényi entropy. For example, one can express the moments of the RDM of|Ψ⟩associated with the regionA = [x_1,x_2]as
    Trłρ^n_A∝̊ ^n⟨Ψ|𝒯(x_1)𝒯̃(x_2)|Ψ⟩^n.
    
    Here|Ψ⟩^n ≡|Ψ⟩^⊗ nis the replicated state, and the non-universal proportionality constant can be absorbed in the normalisation of the twist fields.
    
    
    While Eq. (<ref>) is sufficient to reconstruct the semiclassical predictions of the entropy of particles <cit.> (generated by local operators), when kinks are involved one needs to understand the corresponding generalisation of Eq. (<ref>) to semilocal operators. To the best of our knowledge, this problem has not been considered in the previous literature. We conjecture that the commutation relations betweenμand𝒯read as follows
    𝒯(x)μ^j(y) = μ^j(y)( μ^jμ^j+1·𝒯 )(x),    x>y,
    𝒯(x)μ^j(y) = μ^j+1(y)𝒯(x),    x<y,
    
    with(μ^jμ^j+1·𝒯)(x)the lightest field generated by the fusion
    μ^jμ^j+1×𝒯,
    
    and(μ^jμ^j+1)(x)is a semilocal field of the replica model with disorder lines on replicasjandj+1. Physically, this means thatμcharges the standard twist field𝒯generating a new composite fieldμ^jμ^j+1·𝒯. We give a pictorial representation of Eq. (<ref>) in Fig. <ref>: the branch cut of the twist field is represented as a black dashed line while the disorder line is the red line.
    
    
    
    
    A simple derivation of Eq. (<ref>) can be obtained assuming thatμ(x)is generated by a product of local spin flip operators𝒪(x')inserted at positionx'>xμ(x) = ∏_x'>x𝒪(x).
    
    In the casex>y, a straightforward computation gives
    𝒯(x)μ^j(y) =    𝒯(x)∏_y'>y𝒪^j(y') = ∏_y<y'<x𝒪^j(y') ∏_y'>x𝒪^j+1(y')𝒯(x)=
       ∏_y'>y𝒪^j(y')∏_y'>x𝒪^j(y')∏_y'>x𝒪^j+1(y')𝒯(x) = μ^j(y)μ^j(x)μ^j+1(x)𝒯(x),
    
    where the locality of𝒪(x)has been employed, through Eq. (<ref>), together with the property𝒪(x)𝒪(x)=1. We remark that, while the derivation above can be made rigorous in the lattice, its extension to field theories is not straightforward due to divergences coming from the insertion of fields at the same points; we claim that these issues are only technical and Eq. (<ref>) still holds in the QFT.
    
    In general, the commutation relation between𝒯and the string(μ^j_1…μ^j_k)(x), associated with disorder lines in the replicasj_1…,j_k, can be found in a similar way: each disorder lineμ^jcrosses the branch cut, and it charges the twist field with an additional insertion ofμ^jandμ^j+1. Furthermore, two disorder lines associated with the same replica simplify (due to theℤ_2fusion ruleμ×μ→ 1). In conclusion, via this procedure only twist fields with an even number of disorder lines attached to them are generated, and we refer the interested reader to Ref. <cit.> for the characterisation of their form factors. For our purposes, we only need to recall the properties of the 0-kink form factors, namely the matrix elements between the (2^n) vacua in the replica theory. In particular, due to the topological constraints, it was argued in <cit.> that the only non-vanishing 0-kink form factors for the operators above are
    
    ^n⟨+|𝒯(x)|+⟩^n =  ^n⟨-|𝒯(x)|-⟩^n∝ m^-Δ_𝒯.
    
    Here,m^-1is a correlation length,Δ_𝒯is the scaling dimension of𝒯, and the non-universal proportionality constant depends on the normalisation of the field. A technical observation is in order: these results refer to the theory at infinite size, while here we are interested in the finite-size counterpart. However, as investigated in Ref. <cit.>, we expect deviations from the infinite-size theory to be exponentially suppressed in the system size, and therefore we neglect those in the forthcoming discussion.
    
    
    
     §.§ One-kink entropy
    
    The algebraic relations (<ref>) summarise the essential distinction between (local) particles and kink excitations at the level of entropy of regions. We now claim that these relations are sufficient to give explicit predictions for the Rényi entropies in the semiclassical limit where the regions are large compared to the microscopic lengths. In this section we provide an application of the formalism developed so far by analysing a simple yet non-trivial case, that is the entropy of a single kink in a tripartite geometry. 
    
    We consider an intervalA = [ℓ_1,ℓ_2]and its complementB = B_1∪ B_2 = [0,ℓ_1]∪[ℓ_2,L]and we study the Rényi entropy of the state
    |Ψ⟩ = μ(p)|+⟩.
    
    To do so, we need to compute the ratio
    ^n⟨Ψ|𝒯(ℓ_1)𝒯̃(ℓ_2)|Ψ⟩^n/^n⟨Ψ|Ψ|^⟩n,
    
    where the denominator ensures proper normalisation, and compare it with
    
    ^n⟨+|𝒯(ℓ_1)𝒯̃(ℓ_2)|+⟩^n,
    
    which yields the ground state contribution. Let us first compute the normalisation^n⟨Ψ|Ψ|^⟩nas follows
    
    ^n⟨Ψ|Ψ|^⟩n = ł⟨Ψ|Ψ| ̊⟩^n = ł⟨+|(μ(p))^†μ(p)|+⟩^̊n ∝ L^n,
    
    where Eq. (<ref>) and⟨+|+|=⟩1have been employed.
    
    To proceed with the computation, it is convenient to split the fieldμ^j(p)in terms of its spatial restrictions:
    μ^j(p) = μ^j_A(p)+μ^j_B_1(p)+μ^j_B_2(p).
    
    The operators above have definite commutation relations with the twist fields, which, as a consequence of (<ref>), read
    𝒯(ℓ_1)𝒯̃(ℓ_2)μ^j_B_1(p)    = μ^j_B_1(p)(μ^jμ^j+1·𝒯 )(ℓ_1)( μ^jμ^j+1·𝒯̃ )(ℓ_2),
    𝒯(ℓ_1)𝒯̃(ℓ_2)μ^j_B_2(p)    = μ^j_B_2(p)𝒯(ℓ_1)𝒯̃(ℓ_2),
    𝒯(ℓ_1)𝒯̃(ℓ_2)μ^j_A(p)    = μ^j+1_A(p)𝒯(ℓ_1)( μ^jμ^j+1·𝒯̃)(ℓ_2).
    
    We expand the numerator of Eq. (<ref>) in terms of the restrictions of the fields as
    
    ^n⟨Ψ|𝒯(ℓ_1)𝒯̃(ℓ_2)|Ψ⟩^n    = ^n⟨+| (μ^n(p))^†…(μ^1(p))^†𝒯(ℓ_1)𝒯̃(ℓ_2) μ^1(p)…μ^n(p)|+⟩^n 
       =
    ^n⟨+| (μ^n(p))^†…(μ^1(p))^†𝒯(ℓ_1)𝒯̃(ℓ_2) 
       × (μ^1_A(p)+μ^1_B_1(p)+μ^1_B_2(p))…(μ^n_A(p)+μ^n_B_1(p)+μ^n_B_2(p))|+⟩^n.
    
    We employ Eq. (<ref>), bringing the restrictions ofμ^j(p)to the left of the twist fields, so that they can be eventually contracted withłμ^j(p)^̊†. The contractions produce a total of3^nterms. However, most of these are vanishing in the limit we are considering 
    For instance, since Eq. (<ref>) gives
    
    (μ^j'(p))^†μ^j_C(p) ∝ V_C δ_jj',
    
    for any regionC, it follows that only contractions between disorder lines in the same replica survive. In addition, among the remaining terms, we should only keep those for which the field obtained after performing all the commutation relations is the standard twist field, otherwise, the corresponding vacuum expectation values vanish. 
    In conclusion, only three terms remain from the expansion in Eq. (<ref>), so we can write
    
    ^n⟨Ψ|𝒯(ℓ_1)𝒯̃(ℓ_2)|Ψ⟩^n ≃    ^n⟨+| (μ^n(p))^†…(μ^1(p))^†μ^1_B_1(p)…μ^n_B_1(p) 𝒯(ℓ_1)𝒯̃(ℓ_2) |+⟩^n+
       ^n⟨+| (μ^n(p))^†…(μ^1(p))^†μ^1_B_2(p)…μ^n_B_2(p) 𝒯(ℓ_1)𝒯̃(ℓ_2) |+⟩^n+
       ^n⟨+| (μ^n(p))^†…(μ^1(p))^†μ^2_A(p)…μ^1_A(p) 𝒯(ℓ_1)𝒯̃(ℓ_2) |+⟩^n ∝
        (ℓ_1^n + (ℓ_2-ℓ_1)^n + (L-ℓ_2)^n)×^n⟨+|𝒯(ℓ_1)𝒯̃(ℓ_2)|+⟩^n.
    
    A schematic representation of the mechanism depicted above is sketched in Fig. <ref>.
    
    
    
    
    Putting everything together, we compute the following universal ratio
    ^n⟨Ψ|𝒯(ℓ_1)𝒯̃(ℓ_2)|Ψ⟩^n/^n⟨Ψ|Ψ|^⟩n× ^n⟨+|𝒯(ℓ_1)𝒯̃(ℓ_2)|+⟩^n≃ (r_1^n+r^n_2+r_3^n),
    
    withr_jdefined by
    
    r_1 ≡ℓ_1/L,    r_2 ≡ℓ_2-ℓ_1/L,    r_3 ≡L-ℓ_2/L,
    
    corresponding to the probabilities to find the kink in the regionB_1,A,B_2respectively. We finally express the difference of entropy between|Ψ⟩and|+⟩as
    
    S_n-S_n,0 = 1/1-nlogł r^n_1+r^n_2+r^n_3.̊
    This is the same result (<ref>), up to the vacuum entropyS_n,0that is, in general, non-vanishing in a field theory. Remarkably, we find universality of the entropy difference between the low-lying excited states and the vacuum, which is one of the significant achievements of the field-theoretic framework.
    
    Some comments regarding the regime of small regions are needed. So far, we analysed the scaling limit where the ratios of the sizes are fixed, and they are much larger than the microscopic lengths: in particular, the regimeℓ_2-ℓ_1 ≲ m^-1is not described by the previous formula. Indeed, formally, ifℓ_2=ℓ_1the regionAbecomes the empty set and the entropy vanishes. Unfortunately, this behavior is not recovered by the limitr_2→ 0of Eq. (<ref>), and it is not obvious a priori why there should be a residual entropy when the regionAis much smaller than the rest of the system. The reason is that, albeit the probability that the kink belongs toAvanishes in the limit above, one can still distinguish whether the latter belongs toB_1orB_2just by measuring the local magnetisation inA. In particular, the state becomes locally indistinguishable from the statistical mixture 
    … = r_1 ⟨+|…|+⟩ + r_3 ⟨-|…|-⟩,
    
    and thus one observes an entropy differenceS_n-S_n,0 = 1/1-nlogł r^n_1+r^n_3$̊.
    
    
    
    
    
    
    § KRAMERS-WANNIER DUALITY AND THE ENTANGLEMENT OF ALGEBRAS
    
    
    In the previous sections, we have shown that the entanglement content of kinks explicitly differs from that of particles. This might sound surprising at first, since the ordered/disordered phases are related by duality in the quantum Ising chain: in particular, the kinks in the ordered phase are dual to the particles in the disordered phase. Naively, one could expect that the entanglement content of the corresponding excitations should be the same, but this is not the case. The only reasonable conclusion is that the entanglement of regions is not self-dual under Kramers-Wannier duality. This fact can be understood in terms of the non-local nature of the duality. In particular, it has been shown in <cit.> that the entanglement is transferred from local to non-local degrees of freedom. However, while it is well-established that some local and semilocal operators are related by duality (e.g. the order and disorder operators), the identification of a possible dual of twist fields or of the entropy of a region is less obvious. In this Section, we adopt an algebraic approach, in the spirit of Refs. <cit.>, to answer those questions. The formalism, based on the theory of (finite-dimensional) C^*-algebras, allows dealing systematically with observables that are not local in the spin basis, such as the disorder operators and provides a natural framework to properly define their entanglement content.
    
    We give first a brief review of the Kramers-Wannier duality in the context of the quantum Ising chain. Then, we discuss the notion of density matrix associated with observables in the context of C^*-algebras. Finally, we identify the notion of the duals of twist fields and Rényi entropies associated with regions.
    
    
    
    
     §.§ Kramers-Wannier duality in a quantum spin-1/2 chain
    
    
    The Kramers-Wannier duality, in the context of quantum chains, is associated with the possibility of representing the states of a system using spin or kink variables. Specifically, this proves beneficial for models such as the Ising spin chain, as the Hamiltonian remains local in both representations, albeit ordered phases correspond to disordered phases.
    Here, we are interested in a systematic treatment of the duality map between the two descriptions. While this is usually non-invertible, and the Kramers-Wannier duality is known as a non-invertible symmetry <cit.>, specific sectors associated with given boundary conditions are in one-to-one correspondence as we explain below.
    
    Let us consider the Hilbert space ℋ generated by the 2^L configurations of spins  ± along the z axis, for a chain of length L. This was defined in Eq. (<ref>). We call ℋ the site space, and the configurations of site variables give rise to an orthonormal basis. We consider another Hilbert space ℋ' of dimension 2^L-1, defined similarly:
    ℋ'= Span{|s'_3/2,…, s'_L-1/2⟩,    s'_j+1/2 = ±}.
    
    This is referred to as the bond space,  being generated by configurations of L-1 bond variables. We can associate any site configuration with a bond configuration considering the sign changes of adjacent sites. More precisely, we construct a linear map
    
    T: ℋ→ℋ',
    
    acting on the basis elements as follows
    
    T:|s_1,…, s_L⟩→|s'_3/2,…, s'_L-1/2⟩,    s'_j+1/2 = s_js_j+1.
    
    Roughly speaking, we say that a kink at position j+1/2 is present if there is a change of sign between the sites j and j+1. In the following, we refer to T as the duality map.
    
    We can also consider the adjoint of T as
    
    T^†: ℋ' →ℋ,
    
    such that the associated matrices of T, T^† are hermitian conjugate of each other in the configuration basis. From the definition, it holds T T^† =1, namely T^† is an isometry of ℋ' onto itself. However, T^† T≠ 1; this is not surprising, since the two Hilbert spaces ℋ,ℋ' have different dimensions, and, in particular, T has to display a non-trivial kernel. We fix this issue by focusing on a sector of ℋ such that the associated restriction of T becomes invertible. We choose to fix the last spin and define two sectors of dimension 2^L-1ℋ_±≡Span{| s_1,…, s_L-1,±⟩,    s_j = ±},
    
    satisfying
    ℋ = ℋ_+ ⊕ℋ_-.
    
    From now on we focus on ℋ_+ and, with a slight abuse of notation, we refer to T as the restriction
    
    T:ℋ_+ →ℋ'.
    
    In this way, T becomes invertible and both T^† T and TT^† give the identity operator (on ℋ_+ and ℋ' respectively). 
    
    The duality map we defined above is a correspondence between states,but  we can extend it to a map between observables. For instance, we consider
    End(ℋ_+)    →End(ℋ')
    𝒪   ↦𝒪' ≡ T𝒪T^†.
    
    This is an isomorphism between C^*-algebras, since it is invertible and (𝒪_1𝒪_2)' = 𝒪'_1𝒪_2', (𝒪^†)' = (𝒪')^† hold. One can show, by checking, for example, the action of the operators on basis elements related by duality, that the following relations holds
    σ^z_jσ^z_j+1   ↦σ^z_j+1/2,
    σ^x_j   ↦σ^x_j-1/2σ^x_j+1/2.
    
    Moreover, by representing the one-site operators as products of two-site operators and using the fact that Eq. (<ref>) is an isomorphism, one gets
    σ^z_j   ↦∏_j'≥ jσ^z_j'+1/2,
    ∏_j'≥ jσ^x_j'   ↦σ^x_j+1/2.
    
    This means in particular that local and semilocal operators are mixed with each other by the duality map. It is also worth stressing that, since σ^z_j,σ^x_j generate the algebra End(ℋ_+) (via products and linear combinations), the relations above unambiguously characterise the isomorphism (<ref>).
    
    Finally, we comment on some consequences of the non-locality of the duality map. Let us consider a spatial bipartition ℋ_+ =ℋ_A⊗ℋ_B with
    ℋ_A   = Span{|s_1,…, s_ℓ⟩,    s_j = ±},
    ℋ_B   = Span{|s_ℓ+1,…, s_L-1+⟩,    s_j = ±}.
    
    Here, A and B are two subregions consisting of ℓ and L-ℓ sites respectively. Unfortunately there is no way to define the image of ℋ_A under duality due to the non-locality of the latter: the map T cannot be canonically extended to a map T:ℋ_A→ℋ'_A obtained from a possible spatial bipartition ℋ' = ℋ'_A ⊗ℋ'_B. Nonetheless, it is meaningful to map the observables associated with A to other observables of ℋ'. Namely, the former can be first canonically embedded in the space of the observables of ℋ_+ via
    End(ℋ_A) →End(ℋ_A)⊗1_ℋ_B⊂End(ℋ_+),
    
    and then the duality map (<ref>) can be applied as
    End(ℋ_A)⊗1_ℋ_B→End(ℋ').
    
    The price to pay is that the image of the map above contains observables, as it happens for σ^z_j (with j≤ℓ) in Eq.(<ref>), with a string (product) of σ^z attached to them. In this respect, we say that the isomorphism (<ref>) mixes local and semilocal observables.
    
    
    
     §.§ An algebraic definition of the reduced density matrix
    
    
    Here, we review how to construct the reduced density matrix associated with an algebra. This construction, which is well-established in the context of mathematical physics (see e.g. the textbook <cit.>), is nonetheless mostly overlooked by the vast majority of works on entanglement. For instance, people usually refer to a spatial bipartition of a Hilbert space <cit.>ℋ = ℋ_A ⊗ℋ_B,
    
    and the entanglement of A is probed via the properties of the reduced density matrix ρ_A ∈End(ℋ_A). This approach is satisfactory when the observables of interests are realised locally in the corresponding model. However, in the presence of semilocal operators, or when dealing with maps (such as the Kramers-Wannier duality) that do not preserve locality, a more general approach is needed. 
    
    In the following, we will assume that the algebra of observables 𝒜 is a finite-dimensional C^*-algebra. A state is then a linear functional from 𝒜 to the complex numbers:
    𝒜   →ℂ
    
    a    ↦ a,
    
    and it corresponds to the usual notion of expectation value of the observables in 𝒜. As we show below, it is possible to construct a density matrix ρ_𝒜 as an observable itself requiring
    Trłρ_𝒜 a=̊ a,   ∀ a ∈𝒜.
    
    The trace Trł…$̊ entering the above expression is defined unambiguously thanks to the isomorphism between𝒜and a direct sum of matrix algebras (briefly reviewed in Appendix <ref>), that is a standard result of the representation theory ofC^*-algebras (see e.g. Ref. <cit.>). Specifically, given an orthonormal basis{a_i}of𝒜Trł a^†_i a_j=̊δ_ij,
    
    one expressesρ_𝒜as
    ρ_𝒜 = ∑_j r_j a^†_j,    r_j ≡ a_j.
    
    By applying the definitions above to the bipartition (<ref>) and the algebraEndłℋ_A⊗̊1_ℋ_B⊂End(ℋ), one obtains
    ρ_𝒜 = ρ_A ⊗1_ℋ_B/Trł1_ℋ_B.
    
    The main advantage of this formulation is thatρ_𝒜becomes part of the local algebra𝒜and the local Hilbert spaceℋ_Ais never explicitly used.
    
    In the framework depicted above, it is possible to make sense of the dual of a density matrix under the Kramers-Wannier map considered in the previous Sec. <ref>. Indeed, one can pick the algebra𝒜⊆End(ℋ_+)associated with a region of the site space defined in Eq. (<ref>). Then, given a state…for𝒜and its density matrixρ_𝒜, from the map (<ref>) one can define a dual density matrixρ_𝒜'≡ T ρ_𝒜 T^†⊂End(ℋ').
    
    It is not difficult to show thatρ_𝒜'is just the density matrix of the dual algebra
    𝒜' ≡ T𝒜T^†,
    
    associated with the dual state defined by
     TaT^†' ≡ a,    a ∈𝒜.
    
    Indeed, by making use of the propertyT^† T=1, a straightforward calculation yields
    Trłρ_𝒜' a'=̊Trłρ_𝒜 a=̊ a  =  a'',    a ∈𝒜,
    
    witha' = TaT^†a generic observable in𝒜', that is the defining property of the density matrix. It is worth pointing out that expressions such asTρ_A T^†, withρ_A ∈Endłℋ_A$̊, are meaningless since T is defined at the level of the global Hilbert space only, which is the main reason we adopted the algebraic approach to investigate the duality.
    
    Once the density matrices are defined in the algebraic framework, the usual entanglement measures can be introduced as well. For instance, the Rényi entropy is
    
    S_n = 1/1-nlogTrłρ^n_𝒜,̊
    and it matches the standard definition (in Ref. <cit.>) if applied to the algebra of regions. Moreover, since T is an invertible isometry, from Eq. (<ref>) we get
    Trłρ_𝒜^n =̊Trłρ_𝒜'^n ,̊
    and the two algebras have the same entropy for the corresponding dual states. At this point, it is worth giving a simple paradigmatic example to discuss the consequences of the construction above. We consider a state … of ℋ, defined in Eq. (<ref>), namely the ground state of an Ising chain in its ordered phase. We focus on a region A, we compute the associated entropy, that is the entropy of the algebra End(ℋ_A), then we turn our attention to the Kramers-Wannier duality. For instance, we consider the dual state …', which is a ground state of the paramagnetic phase, and we aim to characterise the entropy of some region in this state. The main problem is that, albeit the states … and …' are related by duality, the algebras of regions are not. In particular, the local observables in the region A are now mapped onto an algebra that is not associated with a region of the system ℋ'. Similar considerations hold for low-lying excited states with a finite number of kinks, dual to states with a finite number of particles.
    
    In conclusion, the Kramers-Wannier duality does not give any direct relation regarding the entropy of regions. Nonetheless, we mention for completeness that, since the duality can be implemented by a Matrix Product Operator (MPO) of bond dimension 2<cit.>, the entropy difference between a state and its dual has an upper bound given by b log 2, with b the number of entangling points.
    
    
    
     §.§ The dual of the twist field
    
    
    In this section, we discuss the twist operators in the algebraic framework and their transformation under Kramers-Wannier duality. In the context of spin chains, these operators have been defined in Ref. <cit.> via their explicit basis representation (see also Ref. <cit.>). Here, as we did for the reduced density matrix, we follow a slightly more abstract construction based on algebraic properties of finite dimensional C^*-algebras, relaxing the hypothesis of the strict locality of the observables in the computational basis. In this way, we can safely discuss the notion of twist operators associated with algebras and their dual.
    
    Let us consider an algebra 𝒜, we replicate it n times and we call 𝒜^⊗ n the replica algebra. 𝒜^⊗ n is generated by the operators a^j, obtained by inserting the elements a ∈𝒜 in the jth replica, for every j:
    
    a^j ≡ 1⊗…⊗ a⊗ 1 ⊗… 1.
    
    The twist operator 𝒯_𝒜∈𝒜^⊗ n is an operator which acts as a replica cyclic permutation j→ j+1. More precisely, one requires
    𝒯_𝒜 a^j = a^j+1𝒯_𝒜.
     
    Eq. (<ref>) is already present in <cit.>, where it is obtained from a more fundamental definition. Instead, here we ask whether this can be truly considered a defining property of the twist field. Clearly, if 𝒯_𝒜 satisfies Eq. (<ref>), then also λ𝒯_𝒜 with λ∈ℂ has the same property, and one may wonder whether (<ref>) fixes the twist operator up to a proportionality constant.
    
    In the discussion below, we assume that the algebra 𝒜 is isomorphic to End(ℂ^d) for a given d; this is the framework of main interest since it corresponds to both the algebra of regions in spin chains, and the (non-local) algebra arising from its Kramers-Wannier dual. In those cases, the only element of 𝒜 commuting with every other element is, up to a constant, the identity 1∈𝒜, and the center of the algebra is trivial. A few minor issues arise for other (finite-dimensional) algebras, and they are discussed in the Appendix <ref>.
    
    The existence of a non-zero operator 𝒯_𝒜 is not obvious a priori (we comment on a general algebraic proof in Appendix <ref>). However, since 𝒜 can be realised as an algebra of matrices, one can employ the explicit construction of Ref. <cit.>, ending up with an operator that has the desired properties. In particular, it is possible to construct 𝒯_𝒜 satisfying
    𝒯_𝒜^⊗ n = Trłρ^n_𝒜,̊
    with …^⊗ n the replica state of 𝒜^⊗ n associated with the density matrix ρ^⊗ n_𝒜∈𝒜^⊗ n.
    
    The uniqueness of the twist operator comes from the fact that 𝒜^⊗ n has a trivial center, since it is isomorphic to End(ℂ^dn), and we give an explanation below.
    We first observe that, from Eq. (<ref>), 𝒯_𝒜𝒯_𝒜^† belongs to the center of 𝒜^⊗ n; therefore, if 𝒯_𝒜≠ 0, it is in an invertible element of 𝒜^⊗ n and its inverse is proportional to 𝒯^†_𝒜. Let us now consider two non-zero twist operators 𝒯_𝒜, 𝒯̂_𝒜 satisfying Eq. (<ref>): (𝒯_𝒜)^-1𝒯̂_𝒜 has to be in the center of 𝒜^⊗ n and, since the latter is trivial, it implies 𝒯_𝒜∝𝒯̂_𝒜. The only arbitrariness in the definition of 𝒯_𝒜 is, therefore, the proportionality constant, which can be fixed unambiguously by requiring the normalisation condition (<ref>).
    
    At this point, we have an operator 𝒯_𝒜 defined canonically at the algebraic level, without explicit reference to any local Hilbert space; in this way, we can give meaning to its dual under Kramers-Wannier duality, as done for the density matrix in the previous section. We start from the isomorphism in Eq. (<ref>), between two algebras. We extend it to an isomorphism of their replica as follows
    𝒜^⊗ n   →𝒜'^⊗ n
    𝒪   ↦𝒪' ≡ T^⊗ n𝒪(T^†)^⊗ n.
    
    Following this procedure, we construct a dual twist operator associated with 𝒜' as
    𝒯_𝒜'≡ T^⊗ n𝒯_𝒜ł T^†^̊⊗ n.
    
    The field above satisfies some natural expected properties. In particular, it allows representing the moments of the dual-density matrix associated with the dual algebra. A direct computation, based on Eqs. (<ref>) and (<ref>), gives indeed
    𝒯_𝒜'^' ⊗ n =  T^⊗ n𝒯_𝒜 (T^⊗ n)^†^' ⊗ n = 𝒯_𝒜^⊗ n = Trłρ^n_𝒜=̊Trłρ^n_𝒜'.̊
    
    As an instructive example, we reinterpret the construction above in the context of the Ising Field Theory of Sec. <ref>. The twist field 𝒯(x) in Eq. (<ref>) is the twist operator associated with the algebra 𝒜 generated by the insertion of ε(y) and σ(y) in the region y>x. In particular, it holds 𝒯(x)σ^j(y) = σ^j+1(y)𝒯(x) if x<y, holds. We can consider another algebra 𝒜', related to 𝒜 via Kramers-Wannier duality, associated with the insertions of ε(y) and μ(y) in y>x. We call 𝒯'(x) the associated (dual) twist operator, which satisfies
    𝒯'(x)μ^j(y) = μ^j(y)𝒯'(x),    y<x.
    
    In other words, 𝒯' seesμ as a local operator. 𝒯 and 𝒯' are explicitly different, and we remark that the commutation relation in Eq. (<ref>) are distinct from those in Eq. (<ref>). In conclusion, there is no reason why the correlation functions of 𝒯, associated with the usual notion of entanglement of the Ising model, should be equal in the ordered and disordered phases: we can only infer relations between 𝒯 and 𝒯'.
    
    
    
    § CONCLUSIONS
    
    
    In this work we studied the entropy of kink excitation, showing that universal results emerge in the limit of large regions: this finding extends previous studies on the entanglement content of quasiparticles. We provided a detailed analysis of specific states of a spin-1/2 chain and made use of the qubit-picture, where computations can be carried out explicitly with elementary methods. Furthermore, we discussed a field-theoretic framework and found that the discrepancy between particles/kinks is ultimately traced back to the algebra of twist fields and the corresponding local/semilocal operators. In particular, the disorder fields  charge non-trivially the twist fields, as shown by Eq. (<ref>), and a new family of twist fields (first introduced and characterised in Ref. <cit.>) is generated. Remarkably, while the qubit picture is simpler to deal with, the field-theoretic approach reveals the universal origin of the entropy difference between the kink states and the vacuum, which is far from obvious.
    
    We pointed out the importance of the notion of entropy of algebras, which is crucial to understanding the lack of correspondence of entanglement entropy under Kramers-Wannier duality. A related, albeit different, discussion found in the previous literature concerns the entropy of the Ising model compared to that of free Majorana fermions: while the two models are related (by Jordan Wigner transformation), a discrepancy in the two-interval entropy has been found in Refs. <cit.>. The origin of the mismatch comes from the computation of the entropy of two distinct algebras: the one generated by local operators (say, Pauli matrices), and the one generated by fermions; this difference can be equivalently expressed in terms of spin-structure in the path integral formalism. We also mention recent results <cit.> on the presence of novel entanglement features in one-dimensional spin models mapped to free fermions, such as universal tripartite information and topological order, previously overlooked.
    
    We remark that the main results of this work do not refer specifically to the Ising model and they do not make use of free fermionic techniques; the only important assumption is the spontaneous symmetry breaking of a ℤ_2 symmetry with kink excitations interpolating between the vacua as a low-lying spectrum. In particular, our field-theoretic approach could be used to discriminate the entanglement content of deconfined kinks and their bound states appearing in the more exotic spectrum of the confining chain (Ising model with both transverse and longitudinal fields) of Ref. <cit.>. Furthermore, we expect generalisations to other systems with distinct symmetry-breaking patterns to be straightforward (e.g. the Potts model in its ordered phase): one should identify the algebraic relations between disordered fields and the twist fields and carry out the contractions to evaluate the correlation functions (as done in Sec. <ref>). In this respect, a promising framework to formulate a general theory, that incorporates both fields with disorder lines and duality lines attached to them (as for Kramers-Wannier) is that of Non-invertible generalised symmetries, recently reviewed in Ref. <cit.>.
    
     Acknowledgements: 
    
    LC acknowledges support from ERC  Starting grant 805252 LoCoMacro. LC thanks Vanja Maric and Michele Fossati for their insightful comments on the manuscript. MM is grateful for funding under the EPSRC Mathematical Sciences Doctoral Training Partnership EP/W524104/1. MM thanks Olalla A. Castro-Alvaredo, David Horvath and Fabio Sailis for the many useful discussions.
    
    
    
    
    
    
    
    
    
    
    § DETAILS OF SOME LATTICE AND QUBIT CALCULATIONS
    
    In this appendix, we collect some technical calculations on the reduced density matrix of multi-kink states on the lattice and within the qubit picture.
    
    
     §.§ RDM of the lattice two-kink state in tripartite geometry
    
    Here, we derive the explicit expression for the reduced density matrix in Section <ref>, that is
    
        ρ^(2)_A_1 ∪ A_3 = 𝒩∑_s_k=±
         k=ℓ_1+1,…,ℓ_1+ℓ_2 ∑_1 ≤ i_1 < j_1 ≤ L-1 
         1 ≤ i_2 < j_2 ≤ L-1   A_2⟨s_ℓ_1 +1,…,s_ℓ_1+ℓ_2|K_+-(i_1)K_-+(j_1)|⟩
        ×   ⟨K_+-(i_2)K_-+(j_2)|s_ℓ_1 +1,…,s_ℓ_1+ℓ_2|_⟩A_2,
    
    with 𝒩≡2/(L-1)(L-2) a normalisation constant. For convenience, we split the indices appearing in the sum into six subsets:
    
        𝒮_1 = {(i,j)|i,j ∈ A_1, i < j }, 
           𝒮_2 = {(i,j)|i ∈ A_1, j ∈ A_2∖{ℓ_1 + ℓ_2}}, 
           𝒮_3 = {(i,j)|i ∈ A_1, j ∈{ℓ_1 + ℓ_2}∪ A_3∖{L}}, 
           𝒮_4 = {(i,j)|i, j ∈ A_2∖{ℓ_1 + ℓ_2}, i < j }, 
           𝒮_5 = {(i,j)|i ∈ A_2∖{ℓ_1 + ℓ_2}, j ∈{ℓ_1 + ℓ_2}∪ A_3∖{L}}, 
           𝒮_6 = {(i,j)|i, j ∈{ℓ_1 + ℓ_2}∪ A_3∖{L}, i < j }.
    
    Then, we compute
    
        A_2⟨s_ℓ_1 +1,…,s_ℓ_1+ℓ_2|K_+-(i_1)K_-+(j_1)|⟩
           = łδ_(i_1, j_1) ∈𝒮_1∏_k ∈ A_2δ_s_k,-|̊%̊s̊⟩̊K_+-(i_1)K_-+(j_1)_A_1⊗|+,…,+⟩_A_3
           + łδ_(i_1, j_1) ∈𝒮_2∏_k = ℓ_1+1^j_1δ_s_k,-∏_k = j_1+1^ℓ_1+ℓ_2δ_s_k,+|̊%̊s̊⟩̊K_+-(i_1)_A_1⊗|+,…,+⟩_A_3
           +łδ_(i_1, j_1) ∈𝒮_3∏_k ∈ A_2δ_s_k,-|̊%̊s̊⟩̊K_+-(i_1)_A_1⊗|K_-+(j_1)⟩_A_3
           +łδ_(i_1, j_1) ∈𝒮_4∏_k = ℓ_1+1^i_1δ_s_k,+∏_k = i_1+1^j_1δ_s_k,-∏_k = j_1+1^ℓ_1 + ℓ_2δ_s_k,+|̊%̊s̊⟩̊+,…,+_A_1⊗|+,…,+⟩_A_3
           +łδ_(i_1, j_1) ∈𝒮_5∏_k = ℓ_1+1^i_1δ_s_k,+∏_k = i_1+1^ℓ_1 +ℓ_2δ_s_k,-|̊%̊s̊⟩̊+,…,+_A_1⊗|K_-+(j_1)⟩_A_3
           +łδ_(i_1, j_1) ∈𝒮_6∏_k ∈ A_2δ_s_k,+|̊%̊s̊⟩̊+,…,+_A_1⊗|K_+-(i_1)K_-+(j_1)⟩_A_3.
    
    In the above expression, to keep the formula as compact as possible we have made the identifications |K_+-(ℓ_1)⟩_A_1≡|+,…,+⟩_A_1, |K_-+(ℓ_1+ℓ_2)⟩_A_3≡|+,…,+⟩_A_3, |K_+-(i)K_-+(ℓ_1)⟩_A_1≡|K_+-(i)⟩_A_1, and  |K_+-(ℓ_1+ℓ_2)K_-+(i)⟩_A_3≡|K_-+(i)⟩_A_3. The RDM is given by the sum of eight blocks, obtained by matching the delta constraints coming from the two matrix elements:
    
        ρ^(2)_A_1 ∪ A_3   = 𝒩∑_(i_1,j_1) ∈𝒮_1
         (i_2,j_2) ∈𝒮_1|K_+-(i_1)K_-+(j_1)⟩_A_1A_1⟨K_+-(i_2)K_-+(j_2)|⊗|+…,+⟩_A_3A_3⟨+,…,+|
           + 𝒩∑_(i_1,j_1) ∈𝒮_6
         (i_2,j_2) ∈𝒮_6|+…,+⟩_A_1A_1⟨+,…,+|⊗|K_+-(i_1)K_-+(j_1)⟩_A_3A_3⟨K_+-(i_2)K_-+(j_2)|
           + 𝒩∑_(i_1,j_1) ∈𝒮_1
         (i_2,j_2) ∈𝒮_6|K_+-(i_1)K_-+(j_1)⟩_A_1A_1⟨+,…,+|⊗|+…,+⟩_A_3A_3⟨K_+-(i_2)K_-+(j_2)|
           + 𝒩∑_(i_1,j_1) ∈𝒮_6
         (i_2,j_2) ∈𝒮_1|+…,+⟩_A_1A_1⟨K_+-(i_2)K_-+(j_2)|⊗|K_+-(i_1)K_-+(j_1)⟩_A_3A_3⟨+,…,+|
           + 𝒩∑_(i_1,j_1) ∈𝒮_3
         (i_2,j_2) ∈𝒮_3|K_+-(i_1)⟩_A_1A_1⟨K_+-(i_2)|⊗|K_-+(j_1)⟩_A_3A_3⟨K_-+(j_2)|
           + 𝒩(ℓ_2 - 1)∑_i_1, i_2 ∈ A_1|K_+-(i_1)⟩_A_1A_1⟨K_+-(i_2)|⊗|+…,+⟩_A_3A_3⟨+,…,+|
           + 𝒩(ℓ_2 - 1)∑_j_1, j_2 ∈{ℓ_1 + ℓ_2}∪ A_3 ∖ L|+…,+⟩_A_1A_1⟨+,…,+|⊗|K_-+(j_1)⟩_A_3A_3⟨K_-+(j_2)|
           + 𝒩(ℓ_1 - 1)(ℓ_2 -2)/2|+…,+⟩_A_1A_1⟨+,…,+|⊗|+…,+⟩_A_3A_3⟨+,…,+|.
    
    
     §.§ Qubit description of three- and four-kink states
    
    
    Here we diagonalise the reduced density matrix of the states with N=3,4 kinks in Eq. (<ref>). These states are
    
        |Ψ^(3)⟩   = √(6 r_1 r_2 r_3)|1^+, 1^-, 1^+⟩
           + √(3 r_1^2 r_2)|2^+,1^+,-⟩ + √(3 r_1 r_2^2)|1^+,2^-,-⟩
           + √(3 r_1^2 r_3)|2^+,+,1^+⟩ + √(3 r_1 r_3^2)|1^+,-,2^-⟩
           + √(3 r_2^2 r_3)|+,2^+,1^+⟩ + √(3 r_2 r_3^2)|+,1^+,2^-⟩
           + r_1^3/2|3^+,-,-⟩ + r_2^3/2|+,3^+,-⟩ + r_3^3/2|+,+,3^+⟩,
    
        |Ψ^(4)⟩   = √(12 r_1^2 r_2 r_3)|2^+, 1^+, 1^-⟩ + √(12 r_1 r_2^2 r_3)|1^+, 2^-, 1^-⟩ + √(12 r_1 r_2 r_3^2)|1^+, 1^-, 2^+⟩
           + √(4 r_1 r_2^3)|1^+,3^-,+⟩ + √(6 r_1^2 r_2^2)|2^+,2^+,+⟩ + √(4 r_1^3 r_2)|3^+,1^-,+⟩
           + √(4 r_1 r_3^3)|1^+,-,3^-⟩ + √(6 r_1^2 r_3^2)|2^+,+,2^+⟩ + √(4 r_1^3 r_3)|3^+,-,1^-⟩
           + √(4 r_2 r_3^3)|+,1^+,3^-⟩ + √(6 r_2^2 r_3^2)|+,2^+,2^+⟩ + √(4 r_2^3 r_3)|+,3^+,1^-⟩
           + r_1^2 |4^+,+,+⟩ + r_2^2 |+,4^+,+⟩ + r_3^2 |+,+,4^+⟩.
    
    We trace out the degrees of freedom of A_2, obtaining the diagonal from of the two RDMs
    
        ρ_A^(3)   = r_2^3 |Φ_0⟩⟨Φ_0| + 3r_1 r_2^2 |Φ_1^-⟩⟨Φ_1^-| + 3r_3 r_2^2 |Φ_1^+⟩⟨Φ_1^+|
           + 6 r_1 r_2 r_3 |Φ_2^-⟩⟨Φ_2^-| + 3r_2(r_1^2 + r_3^2) |Φ_2^+⟩⟨Φ_2^+|
           + r_1 (r_1^2 + 3r_3^2) |Φ_3^-⟩⟨Φ_3^-| +  r_3(3r_1^2 + r_3^2) |Φ_3^+⟩⟨Φ_3^+|,
    
    and
    
        ρ_A^(4)   = r_2^4 |Φ_0⟩⟨Φ_0| + 4r_1 r_2^3 |Φ_1^-⟩⟨Φ_1^-| + 4r_3 r_2^3 |Φ_1^+⟩⟨Φ_1^+|
           + 12 r_1 r_2^2 r_3 |Φ_2^-⟩⟨Φ_2^-| + 6r_2^2(r_1^2 + r_3^2) |Φ_2^+⟩⟨Φ_2^+|
           + 4 r_1 r_2(r_1^2 + 3r_3^2) |Φ_3^-⟩⟨Φ_3^-| + 4 r_2 r_3(3r_1^2 + r_3^2) |Φ_3^+⟩⟨Φ_3^+|
           + 4 r_1 r_3(r_1^2 + r_3^2) |Φ_4^-⟩⟨Φ_4^-| + (r_1^4 + r_3^4 + 6r_1^2 r_3^2) |Φ_4^+⟩⟨Φ_4^+|.
    
    The eigenvectors of ρ_A^(3) are:
    
      * |Φ_0⟩ = |+,-⟩ (no kinks in A)
        
      * |Φ_1^-⟩ = |1^+,-⟩ (one kink in A_1, no kinks in A_3)
        
      * |Φ_1^+⟩ = |+,1^+⟩ (no kinks in A_1, one in A_3)
        
      * |Φ_2^-⟩ = |1^+,1^+⟩ (one kink in A_1, one in A_3)
        
      * |Φ_2^+⟩ = r_1 |2^+,-⟩+r_3|+,2^-⟩/√(r_1^2 +r_3^2) (one magnon in A_1 or one magnon in A_3)
        
      * |Φ_3^-⟩= r_1 |3^+,-⟩+√(3)r_3|1+,2^-⟩/√(r_1^2 +3r_3^2) (3 kinks in A_1 or one kink in A_1 and one magnon in A_3)
        
      * |Φ_3^+⟩= √(3)r_1 |2^+,1^+⟩+r_3|+,3^+⟩/√(3r_1^2 +r_3^2) (1 magnon in A_1 and one kink in A_3 or 3 kinks in A_3),
    
    while those of ρ_A^(4) are:
    
      * |Φ_0⟩ = |+,+⟩ (no kinks in A)
        
      * |Φ_1^-⟩ = |1^+,+⟩ (one kink in A_1)
        
      * |Φ_1^+⟩ = |+,1^-⟩ (no kinks in A_1, one in A_3)
        
      * |Φ_2^-⟩ = |1^+,1^-⟩ (one kink in A_1, one in A_3)
        
      * |Φ_2^+⟩ = r_1 |2^+,+⟩+r_3|+,2^+⟩/√(r_1^2 +r_3^2) (one magnon in A_1 or one magnon in A_3)
        
      * |Φ_3^-⟩= r_1 |3^+,+⟩+√(3)r_3|1+,2^+⟩/√(r_1^2 +3r_3^2) (3 kinks in A_1 or one kink in A_1 and one magnon in A_3)
        
      * |Φ_3^+⟩= √(3)r_1 |2^+,1^+⟩+r_3|+,3^+⟩/√(3r_1^2 +r_3^2) (1 magnon in A_1 and one kink in A_3 or 3 kinks in A_3)
        
      * |Φ_4^-⟩ = r_1 |3^+,1^-⟩+r_3|1^+,3^-⟩/√(r_1^2 +r_3^2) (3 kinks in A_1 and one in A_3 or 1 kink in A_1 and three in A_3 )
        
      * |Φ_4^+⟩=r_1^2 |4^+,+⟩+√(6r_1^2 r_3^2)|2^+,2^+⟩+ r_3^2|+,4^+⟩/√(r_1^4 +r_3^4 + 6r_1^2 r_3^2) (2 magnons in A_1 or 2 magnons in A_3 or 1 magnon in A_1 and one magnon in A_3).
    
    The expressions of the eigenstates show how the superposition mechanism explained at the end of Section <ref> works in detail. It is easy to check that the eigenvalues in Eqs. (<ref>) and (<ref>) match the formulae (<ref>), (<ref>).
    
    
    
    § MORE ON C^*-ALGEBRAS AND TWIST OPERATORS
    
    
    In this appendix, we first review some elementary properties of the finite-dimensional C^*-algebras, and then we discuss some additional details regarding the notions of twist operators and entropy. We remark that infinite-dimensional algebras arise systematically in rigorous treatments of the observables of infinite dimensional spin lattices <cit.> and QFT <cit.>, but their treatment is beyond the purpose of this work. However, in this context, we consider the algebra of finite regions to be finite-dimensional, potentially incorporating appropriate ultraviolet regularisation when dealing with quantum field theories.
    
    The first important result is a classification theorem for those algebras, which goes under the name of Wedderburn–Artin theorem (see for instance <cit.>). It states that every finite-dimensional C^*-algebra 𝒜 is isomorphic to a direct sum of full matrix algebras, namely
    𝒜≃⊕_λEndłℂ^d_λ.̊
    Each term in the previous sum is called a factor, and it is (isomorphic to) an algebra of d_λ× d_λ complex matrices. The center of the algebra 𝒜, defined as
    𝒵(𝒜) ≡{ c ∈𝒜  |  [c,a] = 0   ∀ a ∈𝒜},
    
    is generated by the identity operators of each factor
    𝒵(𝒜) = Span{ 1_λ∈Endłℂ^d_λ}.
    
    In the case of a single factor, which is the usual scenario for the algebra of regions of spin chains, the center is trivial as it is generated by the identity matrix. Nonetheless, in other contexts, as for lattice gauge theories <cit.>, the physical observables may display multiple factors, and we refrain from making specific assumptions in this regard.
    
    An important concept in defining the density matrix and the entropy is the (intrinsic) notion of trace. In the context of algebras, a trace is defined as a positive linear functional
    𝒜   →ℂ,
    
    a    →Trł a,̊
    which satisfies Tr(ab-ba) = 0. Without further hypothesis, it might not be obvious whether a trace exists or it is unique (up to a proportionality constant). The first important result is that, whenever 𝒜 is a full matrix algebra, say
    𝒜 = Endłℂ^d,̊
    then the trace is uniquely defined up to a proportionality constant, and it corresponds to the usual notion of the trace of (d× d complex) matrices. To prove the statement, it is sufficient to observe that the image of the commutator
    
    [a,b]≡ ab-ba,
    
    has dimension d^2-1[
    Given an orthonormal basis {|i⟩} one has [|i⟩⟨j|,|j⟩⟨k|] = |i⟩⟨k| for i ≠ k, and [|i⟩⟨j|,|j⟩⟨i|] = |i⟩⟨i|-|j⟩⟨j|. The linear combinations of those matrices which are obtained via commutators generate the space of traceless matrices.
    ] and it is in direct sum with the span of the identity matrix. Therefore, to define a trace satisfying Eq. (<ref>), it is sufficient to specify its value for the identity operator. The most common normalisation, associated with the usual trace of matrices, is
    Trł 1=̊d.
    
    For algebras with multiple factors, many traces can be constructed in principle, and they are identified by the values of the elements in the center; a natural choice is nonetheless
    Trł 1_λ=̊d_λ.
    
    
    It is worth commenting on the entropy for algebras with multiple factors (see details in Ref. <cit.>). Given a state … for 𝒜 in (<ref>), we split it as
    … = ∑_λ p_λ…_λ.
    
    Here p_λ≡1/ 1_λ is the probability of the factor λ; …_λ is the state conditioned to the factor λ, which satisfies 1_λ_λ =1 and 𝒪_λ =0 if 𝒪∈𝒜_λ' (λ'≠λ). Similarly, the density matrix associated with … is decomposed as
    ρ_𝒜 = ∑_λ p_λρ_𝒜_λ,
    
    with ρ_𝒜_λ the density matrix of …_λ. A simple calculation <cit.> gives the entropy of the state as
    
    S(ρ_𝒜) = ∑_λ p_λ S(ρ_𝒜_λ) - ∑_λ p_λlog p_λ.
    
    We mention that a similar decomposition for the entropy arises whenever the density matrix displays a block-diagonal structure, as pointed out in Refs. <cit.> for U(1) conserving systems. The computation of the moments of the reduced density matrix, and thus of the Rényi entropies, can be given in terms of twist operator for algebras with non-trivial center as well, as we explain below.
    
    First, one defines 𝒯_𝒜_λ as twist operators associated with 𝒜_λ and belonging to 𝒜^⊗ n_λ: no issues arise here since the center of 𝒜_λ is trivial. Then, one introduces
    𝒯_𝒜≡∑_λ𝒯_𝒜_λ∈𝒜^⊗ n,
    
    which acts as a replica cyclic permutation and it does not connect distinct factors. We define the replica state …^⊗ n as the state of 𝒜^⊗ n associated with the density matrix ρ^⊗ n_𝒜∈𝒜^⊗ n. A direct computation gives
    Trł𝒯_𝒜ρ^⊗ n_𝒜=̊   ∑_λ,λ_1,…,λ_n p_λ_1… p_λ_nTrł𝒯_𝒜,λłρ_A,λ_1⊗…⊗ρ_A,λ_n̊̊= ∑_λ p^n_λTrł𝒯_𝒜_λρ^⊗ n_𝒜,λ=
       ∑_λ p^n_λTrłρ^n_𝒜,λ=̊Trłρ^n_𝒜,̊ 
    where we used the fact that contributions associated to distinct factors vanish and only the terms with λ=λ_1=… = λ_n remain.
    
    Finally, we give some remarks regarding the existence of twist operators and their generalisation. In the main text, we mention that a microscopic construction for an operator that implements the replica shift j→ j+1 via Eq. (<ref>) can be provided explicitly. 
    In different contexts, other twist operators are found to be useful in the study of entanglement, as for example composite fields with ℤ_2 lines attached to them, as for the fields in Sec. <ref> and those of Ref. <cit.>. Therefore, a natural question to pose is whether a given transformation of observables can always be implemented via a certain operator. More precisely, given
    ϕ: 𝒜→𝒜,
    
    an algebra automorphism (that is ϕ(ab) = ϕ(a)ϕ(b), ϕ(a^†) = ϕ(a)^†) we look for an element 𝒯^ϕ∈𝒜 satisfying
    𝒯^ϕ a = ϕ(a) 𝒯^ϕ,   ∀ a ∈𝒜.
    
    The cyclic permutation of the replica model in Eq. (<ref>) is indeed a specific example of this general problem. In the case of algebras with trivial central, a general result, under the name of Skolem–Noether theorem<cit.>, ensures that such 𝒯^ϕ always exist: in other words, every automorphism can be realised as an inner automorphism via 𝒯^ϕ a (𝒯^ϕ)^-1 = ϕ(a). This result is far from being trivial, and we are not aware of any simple constructive proof. Nonetheless, a proof of the unicity (up to constants) is definitely simpler, as it follows the same argument of Sec. <ref>, and it relies on the triviality of the center of 𝒜. Lastly, we point out that, in the presence of a non-trivial center, not only the unicity but also the existence of an operator 𝒯^ϕ satisfying (<ref>) is not guaranteed a priori. A simple example is provided by the commutative algebra generated by two central elements a,b (that is, 𝒜≃End(ℂ)⊕End(ℂ)) satisfying a^2=a,b^2=b: the map ϕ(a) = b, ϕ(b) = a is an automorphism, but it clearly cannot be realised as an inner automorphism since the algebra is commutative.
