Resilience of hidden order to symmetry-preserving disorder
Marcello Calvanese Strinati Davide Rossini Rosario Fazio Angelo Russomanno
cond-mat.stat-mech cond-mat.quant-gas quant-ph


 



NEST, Scuola Normale Superiore and Istituto Nanoscienze-CNR, 56126 Pisa, Italy
Département de Physique, Ecole Normale Supérieure / PSL Research University,
CNRS, 24 rue Lhomond, 75005 Paris, France
Department of Physics, Bar-Ilan University, 52900 Ramat-Gan, Israel


Dipartimento di Fisica, Università di Pisa and INFN, Largo Pontecorvo 3, 56127 Pisa, Italy


ICTP, Strada Costiera 11, 34151 Trieste, Italy
NEST, Scuola Normale Superiore and Istituto Nanoscienze-CNR, 56126 Pisa, Italy


NEST, Scuola Normale Superiore and Istituto Nanoscienze-CNR, 56126 Pisa, Italy
ICTP, Strada Costiera 11, 34151 Trieste, Italy




We study the robustness of non-local string order in two paradigmatic disordered spin-chain models, 
a spin-1/2 cluster-Ising and a spin-1 XXZ Heisenberg chain.
In the clean case, they both display a transition from antiferromagnetic to string order.
Applying a disorder which preserves the Hamiltonian symmetries, we find that the transition persists in both models. In the disordered cluster-Ising model we can study the transition analytically – by applying the strongest coupling renormalization group – and numerically – by exploiting integrability to study the antiferromagnetic and string order parameters. We map the model into a quadratic fermion chain, where the transition appears as a change in the number of zero-energy edge modes. 
We also explore its zero-temperature-singularity behavior and find a transition from a non-singular to a singular region, at a point that is different from the one separating non-local and local ordering. The disordered Heisenberg chain can be treated only numerically: by means of MPS-based simulations, we are able to locate the existence of a transition between antiferromagnetic and string-ordered phase, through the study of order parameters. Finally we discuss possible connections of our findings with many body localization. 




Resilience of hidden order to symmetry-preserving disorder
    Angelo Russomanno
    December 30, 2023
==========================================================





§ INTRODUCTION

 

Since Landau, we know that 
phase transitions and symmetry breaking are strictly connected (see for instance Ref. Goldenfeld:book). 
Let us focus on zero-temperature quantum phase transitions (QPTs) <cit.>. 
In the ordered phase, there is a manifold of degenerate ground states which have high entanglement and obey the same symmetries of the Hamiltonian. Nevertheless, any physical ground state breaks this symmetry: because of decoherence, the system always ends up in a superposition of these symmetry-preserving states which has small entanglement and breaks the symmetry. This mechanism manifests in the expectation value Φ of some local operator Φ̂( x) (the order parameter) being different from zero. Symmetry breaking always comes together with long-range order: in the thermodynamic limit, the order parameter in a ground state is infinite-range correlated, 
lim_| x- y|→∞⟨Φ̂^†( x)Φ̂( y) ⟩_ GS=|Φ|^2. 
Although the value of Φ depends on the specific choice of the symmetry-breaking ground state, 
its modulus |Φ|^2 is independent of it. Moreover, this correlator does not depend on the choice of the ground state, even if we consider a non-physical symmetry-preserving ground state. From a physical point of view, the order parameter can be, for instance, the magnetization in the ferromagnetic transition (breaking of rotation symmetry) or the superconducting wave-function in superconductivity (breaking of gauge symmetry). 


This paradigm has been challenged in the last decades by the discovery of topological phase transitions which
are characterized by no local order parameter but by a global rearrangement of the system structure <cit.>.
The landscape of topological phase transitions is extremely rich <cit.>; here we specifically focus on a particular class of topological 
phases occurring in one-dimensional systems,
for which the concept of hidden order has been put forward. In this case, the order parameter 
still exists, but it is a non-local one: it is a string operator involving the system in its globality. 
The most famous example of string order (SO) is in the Haldane phase of a spin-1 isotropic Heisenberg chain <cit.>. 
Let us focus on one-dimensional systems and call 𝒪̂_j,k the SO parameter 
between two sites j and k: long-range order in the thermodynamic limit is given by 
lim_|j-k|→∞⟨𝒪̂_j,k⟩_ GS≠0,
where ⟨ · ⟩_ GS denotes the ground-state (GS) average. 
In this case, an infinite-range non-local correlator is different from zero and does not correspond 
to any local non-vanishing object. 




Examples of transitions to hidden SO phases are the one between a ferromagnet 
and the Haldane phase in spin-1 chains <cit.>, 
and the one occurring in the extended Bose Hubbard model <cit.>. 
In some cases of hidden order the non-local parameter has not yet been recognized or does not exist <cit.>, 
and recently SO has been discovered in a periodically driven spin-1/2 chain <cit.>.

Properties of systems undergoing QPTs are markedly altered by the presence of disorder. 
Disorder shifts the phase transition point and can also create new phases, 
like the Griffiths phase <cit.> in disordered ferromagnets, 
the insulating phase in disordered superconductors <cit.>, the Bose-Glass phase 
in a disordered Bose-Hubbard model <cit.> and the many body localized phase in short-range
interacting quantum systems <cit.>. 
Very remarkable are the works by Dasgupta and Ma <cit.> and by Fisher <cit.>, 
who are able to construct a renormalization group (RG) flow – the “strongest coupling renormalization group” –
to understand phase transitions in such models. This RG method has been used afterwards for a large variety of random quantum and classical systems (see Ref. Igloi_review for a review). Many works focused on disorder and phases 
with local order. The interplay of disorder and non-local SO has been considered in a comparatively small amount of papers. Ref. Igloi focused on disordered spin-1/2 ladders, that may describe properties of spin-1 chains with SO <cit.>. 
The interest on the ladders was due to the fact that the application of strongest coupling RG to spin-1 chains is difficult due to the proliferation of large local spin terms. This difficulty has been overcome in Refs. Boechat1,Boechat,Lajko_PRB; in particular, 


Ref. Boechat found a transition from a Griffiths to a disordered Haldane phase in an S=1 disordered
  antiferromagnetic Heisenberg chain. The numerical demonstration of the persistence of the string order in this model (with a different disorder) and a detailed analysis of the zero-temperature-singularity behaviour has been discussed in Ref. Lajko_PRB. In these works, the form of the disorder has been chosen so to preserve the symmetries of the Hamiltonian. 
 While for small disorder amplitudes the SO is preserved thanks to the topological protection,
  over a given threshold the non-local ordering is broken.
  The topological protection for small disorder is expected from general arguments, however
  the breaking of SO for strong disorder and the properties of this transition are not trivial at all <cit.>.



As far as we know, there are no studies on the effect of a strong disorder on the properties of a transition from a string non-local order
to a local order (like ferro- or antiferro- magnetic).
Here, we make a progress in this direction considering disordered spin-1 and spin-1/2 chains 
whose clean counterparts (the XXZ Heisenberg chain and the cluster-Ising model) are well known to display 
a transition from an antiferromagnetic (AF) to a string-ordered phase (the isotropic Heisenberg chain 
is the Haldane model). 
In the first case we resort to numerics: using a variational algorithm based on the formalism
of matrix-product states (MPS) <cit.> we find that a phase with SO actually does exist, 
even in the disordered model. 
We locate the phase transition point and observe that it is shifted by the disorder, with respect to the clean case. 
In the case of the disordered cluster-Ising model we can go further: by applying 
a strongest coupling RG transformation, we show the existence of an AF and a string-ordered phase, 
and analytically find the phase transition point. 
In this case, the model is solvable because it can be mapped to a disordered non-interacting fermion chain. 
In the fermionic representation, the transition is of topological nature and is 
characterized by the appearance of zero-energy edge modes.
Finally, we 
  study the zero-temperature-singularity behavior of the model, by numerically analyzing the properties
  of the logarithmic gap distribution and of the inverse average of the gap <cit.>.
  We find a transition between a non-singular and a singular behavior when the disorder strength is increased.
  This transition point is different from the one between SO and AF: for our specific form of disorder, we see the existence
  of a non-singular and a singular SO phase, similarly to the findings of Ref. Lajko_PRB.
  Moreover we observe that the AF phase is singular.



The paper is organized as follows. In Sec. <ref>, we introduce the two models we are considering. 
In Sec. <ref>, we study the AF/SO transition in the cluster-Ising chain. 
We do this both analytically (through the strongest coupling RG) and numerically (evaluating the appropriate 
correlators thanks to the free-fermion mapping); we also study the properties of edge modes 
in the fermionic representation. 
In Sec. <ref>, we study the AF/SO transition in the Heisenberg XXZ chain, 
by resorting to MPS numerical simulations. 
Finally, in Sec. <ref>, we draw our conclusions. 
In Appendix <ref>, we show details on the strongest coupling RG transformation applied 
to the cluster-Ising chain, and in Appendix <ref> we discuss the appearance 
of kinks in the MPS approximation of the GS, which are due to a numerical artifact.




§ THE MODELS



In the following, we will consider two spin chains that exhibit a zero-temperature QPT
between a Néel-like AF phase and a phase that is characterized by non-local SO. They are the spin-1/2 cluster-Ising model
and the spin-1 XXZ Heisenberg chain. The two models enjoy a local ℤ_2 and 𝔻_2 symmetry, respectively. We thus expect the SO phase to exist also in the presence of disorder, as long as the local symmetry is not broken <cit.>. In this section, we introduce the two models by presenting their Hamiltonian together with some of their properties and the relevant antiferromagnetic and string order parameters.





 §.§ Spin-1/2 cluster-Ising model
 

The simplest and exactly solvable model in this context is the so called cluster-Ising model (CIM) <cit.>,
which is described by the Hamiltonian


    Ĥ_ CIM = ∑_j [ - J_j σ̂_j-1^xσ̂_j^zσ̂_j+1^x
      + λ_j σ̂_j^yσ̂_j+1^y ]  .


Here, σ̂^α_j (with α = x,y,z) denote the spin-1/2 Pauli matrices on the j-th site of the chain
(j=1, …, L, where L is the chain length),
while J_j and λ_j denote the (possibly site-dependent) three-spin and two-spin coupling terms, respectively.
In the following, we will always adopt open boundary conditions (OBC).
After a standard Jordan-Wigner transformation <cit.> of the form
ĉ_j = ( ∏_m=1^j-1σ̂^z_m ) σ̂^-_j/2,
the Hamiltonian in Eq. (<ref>) can be mapped onto a free-fermion model with nearest-
and next-to-nearest-neighbor hopping terms


    Ĥ_ CIM = ∑_j [     - J_j (j-1-j-1)(j+1+j+1)
       + λ_j (j+j)(j+1-j+1) ]  ,


in a way similar to the one-dimensional quantum Ising model in transverse field <cit.>. 
This Hamiltonian is integrable, and its dynamics can be reduced to that of non-interacting 
fermionic quasiparticles: the GS is the BCS state without quasiparticles.

It has been proven <cit.> that, in the homogeneous case (i.e. J_j=J and λ_j = λ for all j),
the system features a QPT between a conventional AF phase along the y direction and a phase with 
non-local SO (the so called “cluster phase”), when decreasing λ/J and crossing the critical value of 1.

The y-antiferromagnet is detected by a local order parameter. In agreement with the discussion 
in the Introduction, we can express it as an infinite-range correlator:


    ( 𝒮_[1/2]^y)^2 = 1/4 lim_l→∞ (-1)^l  ⟨σ̂_k^y σ̂_k+l^y ⟩ .


To be precise, here we are considering the square modulus of the order parameter, 
which is the same for all the degenerate symmetry-breaking ground states. 
Since this model breaks the ℤ_2 symmetry (see below for more details), there are exactly 
two symmetry-breaking ground states which differ for the sign of the order parameter.
Moreover, because we are evaluating the modulus of the order parameter as the limit of a correlator, 
it is not important to select the symmetry breaking ground states: the correlators are the same 
on all the states of the GS manifold. From a technical point of view, the mapping 
to the fermionic model is very important in order to evaluate this correlator. 
Thanks to Wick's theorem, it can be expressed as a Toeplitz determinant <cit.> 
of specific single-particle fermionic correlators 
of the Hamiltonian (<ref>) – see Subsec. <ref> for more details. 

On the opposite, the cluster phase is characterized by a non-vanishing value of the non-local SO parameter, which is expressed as the infinite-range limit of a non-local correlator:


    𝒪_[1/2]^z = lim_l →∞ (-1)^l  ⟨σ̂_k^xσ̂_k+1^y 
      [ ∏_n=k+2^k+l-2σ̂_n^z ] σ̂_k+l-1^y σ̂_k+l^x ⟩ .


From a technical point of view, this non-local correlator is evaluated by applying to the Hamiltonian 
a duality transformation (see Ref. Smacchia_2011 and Appendix <ref> for details). 
This non-local unitary transformation maps the Hamiltonian onto itself, with λ_j exchanged with J_j 
[see Eq. (<ref>)]. Moreover, the correlator of Eq. (<ref>) 
is mapped on the AF correlator of Eq. (<ref>), which can be evaluated as explained above.

Unless specified, here and in the following equations, we will omit the subscript ⟨ … ⟩_ GS 
and consider all the expectation values on one of the degenerate ground states of the system. 
As remarked before, it is not important which GS in particular, 
being the correlators independent of the specific choice of the state in the GS manifold. 
Moreover, in the presence of disorder, we will average over many realizations: 
this operation will be denoted by an overline …. 
We notice that the disorder is not translationally invariant: 
translation invariance is restored after averaging over the disorder realizations.

From a numerical point of view, we can only evaluate finite range correlators: 
𝒪_[1/2],l^α (with α=x,y,z) denotes the string correlator over a range l,
such that 𝒪_[1/2]^α = lim_l →∞𝒪_[1/2],l^α.
In a similar way, we define the finite-range AF correlator as
( 𝒮_[1/2],l^α)^2 such that
( 𝒮_[1/2]^α)^2 = lim_l →∞( 𝒮_[1/2],l^α)^2. 

Let us now focus on the version of the Hamiltonian with disorder, which we are studying in the rest of the work.
A necessary condition to preserve the string-ordered phase is the choice of a disorder that preserves 
the symmetry of the model: this is the symmetry broken by the GS 
in the phase transition. As mentioned above, the Hamiltonian in Eq. (<ref>) enjoys a ℤ_2 symmetry, 
since it is invariant under a rotation of angle π around the z axis:
V̂^†Ĥ_ CIMV̂ = Ĥ_ CIM with 
V̂ = exp(-iπ/2∑_jσ̂_j^z).
Being the symmetry group generated by only two operators (V̂ and 𝕀̂), the GS manifold
can have at most dimension two, so there can be at most two symmetry-breaking ground states.
The invariance under the operator V̂ implies:


    [ ∏_jσ̂_j^z,Ĥ_ CIM]=0 .


A possibility to satisfy this condition in a disordered situation is to assume that both the site-dependent two-spin and
the three-spin couplings are taken randomly from some probability distribution. In the following, we will specifically address the situation in which
both the λ_j and the J_j are uniformly distributed over some interval
λ_j ∈ [0,λ_ max] and J_j ∈ [0, J_ max], for all j=1, …, L.





 §.§ Spin-1 XXZ Heisenberg chain
 

The other model that we are going to focus on is a spin-1 XXZ Heisenberg chain, given by the Hamiltonian


    Ĥ_ XXZ = J ∑_j [ Ŝ^x_jŜ^x_j+1+Ŝ^y_jŜ^y_j+1+Λ_jŜ^z_jŜ^z_j+1]   ,


where Ŝ_j^α (with α = x,y,z) now denotes the spin-1 operators on the j-th site, J is the energy scale
of the nearest-neighbor spin coupling, and Λ_j is the anisotropy factor along the z axis, at site j.
Contrary to the CIM, the spin-1 Heisenberg chain is a non-integrable model and cannot be easily diagonalized.
For this reason, numerical approaches based on exact diagonalization techniques or on MPS
are usually employed in order to capture the GS physics.

In the clean case, that is for Λ_j = Λ,  ∀ j, model (<ref>) is known to display 
a phase transition between a topological phase, usually referred to as the Haldane phase <cit.>, 
and a non-topological Néel AF phase.
Such transition has been studied in some details in the literature, and is expected to occur
for Λ≈ 1.186… <cit.>.

In an open chain, the GS in the Haldane phase is identified by the presence of gapless spin-1/2 modes on top of a gapped bulk, which make the GS
four-fold degenerate in the thermodynamic limit. This phenomenology is related to a hidden 𝔻_2 symmetry breaking, described by the dihedral group of rotations <cit.>


    𝒢_𝔻_2 = {Î,e^iπ∑_nŜ^x_n,e^iπ∑_nŜ^y_n,e^iπ∑_nŜ^z_n}  .


Because of the presence of such edge modes, a state in the Haldane phase is characterized by a hidden long-range order <cit.>,
which cannot be revealed by any expectation value of simple two-point correlators
⟨Ŝ^α_kŜ^α_k+l⟩.
Indeed, their GS expectation values vanish in the limit l →∞.
This hidden order can be seen by defining a non-local SO parameter for a spin-1 chain 
in a way similar to what has been done for the CIM (<ref>):


    𝒪^α_[1] = lim_l →∞𝒪^α_[1],l =  lim_l →∞⟨Ŝ^α_k [ ∏_n=k+1^k+l-1e^iπŜ^α_n ] Ŝ^α_k+l⟩ ,


such that the system is said to posses SO if the above limit is finite and nonvanishing.

The presence of hidden order can be understood by remapping the model onto a ferromagnetic chain 
with four symmetry-broken states,
by means of the Kennedy-Tasaki transformation <cit.>: string correlators are mapped
onto two-point correlators, which indeed reveal the presence of ferromagnetic order. 
This is similar to what happens in the CIM, where a duality transformation maps the SO correlator 
of Eq. (<ref>) onto the AF correlator of Eq. (<ref>). 

The physical meaning of the SO parameter (<ref>) is very clearly expressed in Ref. Scalapino; 
we briefly review it here for reader's convenience: if we measure the value of the spin projections along α, 
we can find 1,0 or -1, being this a spin-1 chain. Since the GS is a superposition of spin eigenstates, 
each time we measure we find a different sequence of +1,0 or -1, with some probability. 
The fact that the operator (<ref>) has a non-vanishing expectation value means that, if we withdraw the zeros 
from any of these sequences, we find alternatively +1 and -1: without the zeros, the system behaves 
antiferromagnetically. This property cannot be witnessed by any local operator: only the non-local string operator can do.

On the other hand, the presence of long-range AF order along the z axis in the Néel phase is witnessed 
by a nonzero value of the staggered two-point correlator 


    (𝒮_[1]^α)^2 = lim_l →∞( 𝒮_[1],l^α)^2 
      = lim_l →∞ (-1)^l  ⟨Ŝ^α_k Ŝ^α_k+l⟩ ,


with α = z (this is the square value of the AF order parameter along α).
Thus, we can identify the Haldane phase of the XXZ chain by a nonzero SO [Eq. (<ref>), for all α]
and a vanishing expectation value of the staggered correlator in Eq. (<ref>), for any α.
Conversely, the Néel phase is identified by a vanishing SO 𝒪^α_[1] for α=x,y,
and by a nonzero value of the staggered AF order parameter 𝒮^z_[1].
We point out that we cannot use the SO along z as an order parameter of the Haldane phase,
because the observable 𝒪^z_[1] is nonzero both in the Haldane and in the Néel phase <cit.>.
This can be simply understood in the large-Λ limit, where the GS is given by a product
of consecutive states with opposite spin projection. 
From Eq. (<ref>), it is clear that 𝒪̂^z_[1],l evaluated on such GS is exactly -1
for all l, and in particular for l →∞.

We will model the disorder by taking Λ_j as a random variable,
which is uniformly distributed between Λ_ min and Λ_ max.
We stress that a necessary condition for the GS to possess SO is to enjoy a unitary and local symmetry <cit.>. 
As stated before, since the presence of non-uniform anisotropy Λ_j
does not break the 𝔻_2 symmetry, the SO phase is expected to be present also in the disordered
XXZ model, at least for some range of values Λ_ min and Λ_ max (see Appendix <ref> for an example of destruction of string order when a symmetry breaking term is added to the Hamiltonian).





§ ANTIFERROMAGNETIC – CLUSTER PHASE TRANSITION IN THE DISORDERED CIM
 


 §.§ Strongest coupling RG approach
 

To understand whether the CIM can undergo a QPT or not, we perform a strongest coupling RG analysis 
in the thermodynamic limit L →∞, very similar in spirit to the one used in Ref. Fisher_1995 
for the disordered transverse field Ising model. We consider the largest value of the coupling


    Ω_I=max{J_j, λ_j} .



The idea is to diagonalize the part of the Hamiltonian related to this coupling, 
assuming that it is so large that the corresponding term of the Hamiltonian can be considered as non-interacting, in a crudest approximation. The rest of the chain
can be considered as a perturbation which changes the GS energy of this subsystem at second order in the coupling. 
At the end of the renormalization step, we take only the perturbed GS of the renormalized subsystem, 
and discard the rest of its local Hilbert space. In this way, at each renormalization step, we reduce 
the energy scale at which we are looking at the system: the considered site is renormalized away 
and an effective low-energy coupling term is generated. After many applications of the renormalization step, 
we asymptotically reach the GS and we can find its properties.
One can distinguish between two possible cases, depending whether the largest coupling is a given J_j, 
or a given λ_j. We provide the derivation in full detail in Appendix <ref>; 
here, we focus on its main points and their physical meaning.

We start assuming that the largest coupling is J_j: the part of the Hamiltonian corresponding to it is given by


    Ĥ_0 = -J_jσ̂_j-1^xσ̂_j^zσ̂_j+1^x  ,

while the coupling to the rest of the system can be described by the following operator:


    V̂=λ_j-1σ̂_j-2^yσ̂_j-1^y+λ_j+1σ̂_j+1^yσ̂_j+2^y  .


As detailed in Appendix <ref>, it is possible to treat the term V̂ perturbatively, 
applying a first-order perturbation theory to the four-fold degenerate GS of Ĥ_0.
After diagonalizing the resulting second-order perturbation matrix,

we project over one of the perturbed ground states, ending up into eliminating site j and generating a new coupling


    -λ_jσ̂_j-2^yσ̂_j+2^y   with  λ_j≃λ_j-1λ_j+1/J_j .


The operators σ̂_j-2/j+2^y in Eq. (<ref>) are in principle different from the unrenormalized ones:
they coincide with them up to quartic terms in λ/J.

In the opposite case, where the largest coupling is λ_j, one can apply to Eq. (<ref>) 
a duality transformation


    μ̂_j^x = ∏_k=1^jσ̂_k^z  ,     μ̂_j^z = σ̂_j^xσ̂_j+1^x ,


which maps the Pauli operators σ̂_j^α onto different Pauli operators μ̂_j^α. 
After the application of this transformation, the CIM Hamiltonian in terms of σ̂_j^α [Eq. (<ref>)] is re-expressed in terms of μ̂_j^α [see Eq. (<ref>) in Appendix <ref>]: the transformed Hamiltonian has the same form of Eq. (<ref>), but λ_j and J_j are now exchanged.
The term with the largest coupling which has to be renormalized is now
Ĥ̃̂_0 = -λ_jμ̂_j-1^xμ̂_j^zμ̂_j+1^x.
Applying to it the same analysis of the first case, we see that the renormalization procedure eliminates the site j
in the dual representation, and generates the term


    -J_j+1μ̂_j-2^yμ̂_j+2^y   with  J_j=J_j-1J_j+1λ_j-1 .



In the limit of many RG steps, and after applying the central limit theorem, it is possible to see that 
the disordered CIM is equivalent to a system with couplings 


    λ    =    exp[2l (  logλ - log J ) ] 
    J    =    exp[ 2l (  log J - logλ) ]   ,


where l is the number of consecutive 
renormalized sites, which in principle can be different for different sites of the renormalized model (all the details of the calculation are in Appendix <ref>). 
We can distinguish three cases: 

I: logλ>log J. In this case, λ 
is larger than J exponentially in the number of iterations of the renormalization step. 
In the limit of infinite iterations, J is vanishingly small with respect to λ: 
only the AF terms survive. Therefore, the RG flows to an AF condition and the system is antiferromagnetic. 


II: logλ<log J. Similarly to case I, 
in the limit of infinite iterations, λ is vanishingly small with respect 
to J: only the three-body (cluster) terms survive. 
Looking at the problem in the dual representation, we see that the dual system has only the AF term 
and then the RG in this representation flows to an antiferromagnetic condition. 
Going back to the original representation, we see that the system flows to the SO phase. 


III: logλ=log J. Here, the RG flows to a uniform system 
with λ=J = 1: the low-energy behavior of the model is equivalent 
to a uniform model at the critical point between the SO and the AF phase. 

Therefore, also in the disordered model, we can predict a transition between AF and SO phase,
occurring for 


    logλ = log J .


For logλ>log J the model is AF, while for 
logλ<log J it displays SO.
These results are very similar to those found in the transverse field Ising model <cit.>.





 §.§ Numerical analysis of the two phases
 

Just to fix the ideas, let us now analyze the case in which there is no disorder on J (J_j = 1,  ∀ j),
and each λ_j is taken from a uniform distribution among 0 and some λ_ max.
In this case, the condition in Eq. (<ref>) implies that the transition point is located at


    λ_ max,c^(∞) =  e≈ 2.718… ,


e being the Neper number. The superscript “(∞)” in Eq. (<ref>) denotes that this is the critical point in the limit L→∞.
By exploiting the fact that the CIM is exactly solvable, we can explore the behavior of the long-range string and 
AF correlator for considerably long system sizes, after averaging over an ensemble of several disorder realizations. 

Here, we recall that, fixing the realization of the disorder, a Jordan-Wigner transformation is able to map 
Eq. (<ref>) into a free-fermion diagonal form (irrespective of the presence or absence of translation invariance):


    Ĥ=∑_μ=1^Lϵ_μ(2 γ̂^†_μγ̂_μ - 1 ) ,


where the single quasi-particle operators γ̂_μ are defined in terms of the local fermionic
operators ĉ_j according to


    γ̂_μ = ∑_j=1^L(U_jμ^ *ĉ_j + V_jμ^ *ĉ^†_j ) ,


and U_jμ^ *, V_jμ^ * are the coefficients of the 2L× 2L unitary matrix which
diagonalizes the appropriate 2L× 2L Hermitian matrices forming the Hamiltonian 
(see for instance Ref. Russomanno_JSTAT13 for more details on this method). 
The GS is the one which is annihilated by all the γ̂_μ operators (it has a BCS form). 
Thanks to this property, we can evaluate the AF correlator in Eq. (<ref>).
Applying Wick's theorem to the BCS Gaussian state, we can write the AF correlator as a Toeplitz determinant:


    (𝒮_[1/2],l^α)^2=
      1/4 |[     G_k,k+1           ⋯     G_k,k+l;           ⋮                       ⋮; G_k+l-1,k+l           ⋯ G_k+l-1,k+l ]|  ,


where, for each disorder realization, we have defined the two point fermionic correlators
on the GS corresponding to that realization:


    G_j,m = ⟨ (j-j)(m+m) ⟩ .


Inverting Eq. (<ref>) and using the fact that the GS is annihilated by all the γ̂_μ, 
we can evaluate this correlator as


    G_j,m=∑_μ(V_jμ^ *-U_jμ)(U_mμ^ *+V_mμ) .


The SO parameter in Eq. (<ref>) is evaluated by applying the duality transformation (<ref>), 
which maps it onto an AF correlator of the form in Eq. (<ref>), and the Hamiltonian onto 
another Hamiltonian of the same form. Finally, the results obtained through these formulas 
are averaged over N_ av realizations of disorder.

The outcomes of our computations for a given finite size are reported in Fig <ref>. 
On the upper panel we plot 𝒪^z_[1/2],l and (𝒮^y_[1/2],l)^2. 
In order to avoid unwanted boundary effects, we evaluate the correlators between sites 
that are far away from the chain ends (see the caption for details).
These quantities approximate the order parameters of Eqs. (<ref>) and  (<ref>), 
since in numerical simulations we can consider large, but yet finite system sizes. 
We see in the upper panel of Fig <ref> that, when 𝒪^z_[1/2],l vanishes, 
(𝒮^y_[1/2],l)^2 appears with a crossover: 
this is an indication that there is a transition from a string-ordered to a y-antiferromagnetic phase 
in the thermodynamic limit, even in the presence of disorder. 





We can also analyze fluctuations over the disorder of the two order parameters. 
As before, we calculate finite-range correlators. Focusing on the AF order, 
we can define the fluctuation as


    δ (𝒮^y_[1/2],l)^2 = [ ⟨σ̂_k^y σ̂_k+l^y ⟩^2-⟨σ̂_k^y σ̂_k+l^y ⟩^2 ]^1/2  ,


where the expectation value has to be intended over the GS of any specific disorder realization. 
The definition for the SO fluctuation (δ𝒪^z_[1/2],l) is analogous, 
after replacing the correlator of Eq. (<ref>) with that of Eq. (<ref>). 
The results are shown in the lower panel of Fig. <ref>. We can see that fluctuation is different from zero 
only when the corresponding order parameter is nonvanishing (compare with the upper panel): 
the finite-size crossover appears also in the behavior of fluctuations.


Until now, we have considered signatures of the transition in finite-length correlators. 
As we can observe in Fig. <ref>, finite-size effects are evident: for a given system size, 
we actually see a crossover, and there is a region where both the finite-range order parameters are different from zero. 
Moreover, if we identify the transition with the point where the curves of the two finite-range order parameters cross, 
we get a result that is different from the theoretical prediction (<ref>). 
In order to properly infer the behavior in the thermodynamic limit L,l→∞, 
we perform a finite-size scaling analysis. 

In order to reduce the effect of the fluctuations induced by the noise, we need to perform a coarse graininig in 1/l of disorder-averaged correlators. 
More precisely, we proceed in the following way. We fix the value of L and, in order to avoid finite-size boundary effects, we fix an appropriate l_0 and consider the disorder-averaged correlator between the site l_0 and the site l with l_0+l varying between 0 and L - 2l_0. Then, we coarse grain this correlator: we consider the interval in which the quantity 1/l varies, divide this interval in windows of width δ(1/l) and perform the average of the correlator over each window. For each of the resulting values we evaluate the uncertainty as 
the maximum over the corresponding window of the disorder fluctuation of the correlator: applying the central limit theorem, this fluctuation is given by the value in Eq. (<ref>)
divided by √(N_ av). 

We label each of the windows over which we coarse grain with its central value 1/l: for each value of 1/l we locate the approximate transition point  as the value of λ_ max where the two coarse-grained correlators cross; we denote the crossing point as λ_ max,c^(l). In the upper panel of Fig. <ref> we show λ_ max,c^(l) versus 1/l: taking into account the error bars, we see a behaviour consistent with a convergence towards

the theoretical value of Eq. (<ref>) when 1/l→0. 
The error bars are evaluated in the following way: each time we add or subtract the uncertainties discussed above to the coarse-grained disorder-averaged 

correlators and then we  take the half-dispersion of the 4 resulting estimates of the crossing. 



The lower panel of Fig. <ref> shows that the height w^(l) of the crossing point tends to zero when 1/l→0 (the error bars are evaluated with the same method used for λ_ max,c^(l)). 
We have checked this fact by fitting with a straight line a·(1/l), and found 
 (a=317± 4)×10^-5 (red line in Fig. <ref>).
This means that, in the thermodynamic limit, w^(∞)=0, and thus there is no region where both the order parameters 
are non-vanishing: when one vanishes the other appears, as appropriate for a QPT. 










 §.§ Edge modes


For a clean system, the QPT in the spin chain maps to a topological transition in the fermionic representation. 
In the case of periodic boundary conditions, the AF phase corresponds to winding number one in the fermionic picture, 
while the SO phase corresponds to winding number two <cit.>. <cit.>
Taking OBC, the topological nature of the system appears through the existence 
of zero-energy boundary modes <cit.>: diagonalizing the fermionic Hamiltonian, 
some vanishing ϵ_μ appear in Eq. (<ref>). 
For each phase, there is a fixed number of zero-energy modes, and their amplitudes U_j and V_j 
[see Eq. (<ref>)] are localized on the edges of the system. 
The AF phase displays one zero-energy mode, while the SO phase has two zero-energy modes (edge modes
in uniform fermionic Hamiltonians very similar to Eq. (<ref>) have been studied in Refs. prb85,prb88).

Even in the presence of disorder, we numerically observe the persistence in the spectrum of zero-energy modes
(two modes in the SO phase, and one mode in the AF phase). 
Two examples of this fact are reported in Fig. <ref>. 
Here, we choose a specific disorder realization and show the single quasi-particle spectrum ϵ_μ 
for a case where the system shows SO (λ_ max=0.8) and a case where it is AF (λ_ max=4.8). 
In the first situation there are two levels with energy many orders of magnitude smaller than the others; 
in the second one there is a single level with this property: these levels correspond to the boundary modes discussed above 
(the energy is not exactly zero, due to the numerical round-off errors). 
We have verified that the same structure of the spectrum appears for any disorder realization. 

These edge modes are topologically protected, since they only depend on global properties of the system: 
they cannot be destroyed by local perturbations (like disorder) if the perturbation is weak enough. 
That is why, if we add disorder, two edge modes and the associated SO persists for λ_ max small, 
and one edge modes and AF order persist for λ_ max large. 
For λ_ max around the transition, the disorder is strong enough to move the transition point.







 §.§ Thermodynamic singularities


Disordered systems can display phases where the thermodynamic quantities show singularities in the limit of vanishing temperature <cit.>. This can be seen from the behaviour of the distribution of the energy gap Δ of the Hamiltonian.
If for small Δ the logarithmic energy gap distribution behaves as P(logΔ)∼Δ^1/z (z is the so-called dynamic exponent), it
is easy to show that the excitation energy over the ground state behaves as E_ ex(T)∼ T^1+1/z at low temperatures. Therefore, the low-temperature specific heat
behaves as C∼ T^1/z and its derivative shows a divergence in the limit of T→ 0 when z>1. The ranges of parameters where this happens are called singular regions <cit.>; in order to find them we numerically consider the properties of P(logΔ) and we check that it behaves as a power law for small Δ (see some examples in Fig. <ref>).



Applying a linear fit to the plots of log[P(logΔ)] vs logΔ we are able to estimate the value of the dynamic exponent z which we plot in Fig. <ref>. We see that there is a singular region with z>1 for λ_ max above a threshold λ_ max^c s≃2. In order to have a better estimate of λ_ max^c s, we follow Ref. Lajko_PRB and consider the behaviour of the inverse average of the gap defined as


    Δ^ iv≡[(1/Δ)]^-1 ,


where the average is performed over the disorder distribution. This object vanishes whenever the system is in a singular region <cit.> with z>1; we show results for our case in Fig. <ref>. We see that Δ^ iv vanishes for λ_ max>λ_ max^c s=2, confirming that in this parameter range our system is singular. Moreover, we can numerically find that Δ^ iv vanishes as a power law when λ approaches the transition point λ_ max^c s from below: we have Δ^ iv∼(λ_ max^c s-λ_ max)^μ_Δ with μ_Δ=2.06±0.01.

In conclusion, we see a transition to a singular regime which occurs at a critical value λ_ max^c s different from the critical point λ_max, c separating the SO and the AF phase. While for λ_ max>λ_ max, c the system is AF and singular, we have a non-singular  SO phase (z<1) for λ_ max<λ_ max^c s and  there is a singular SO phase (z>1) for λ_ max^c s<λ_ max<λ_max, c. This behaviour is strictly reminiscent the disordered S=1 antiferromagnetic Heisenberg chain <cit.> where there is a gapped Haldane phase (z<1) and a singular Haldane phase (z>1), both showing SO. We emphasize that the singularity structure of the phases
is strictly related to the specific form of the disorder: for instance, taking λ uniform and J_j uniformly distributed between 0 and some J_ max, we would have seen a non-singular AF phase together with a singular AF phase, while the SO phase would have been fully non-singular.








§ HALDANE TO NÉEL PHASE TRANSITION IN THE DISORDERED SPIN-1 XXZ MODEL
 

We now switch to study the zero-temperature properties of the disordered spin-1 XXZ Hamiltonian 
of Eq. (<ref>). 
Since this model is non integrable, in order to find the GS of a given finite-size system,
we resort to a variational search on the class of MPS <cit.>.
In our simulations, we analyze chains of up to L=240 sites and
choose a maximum bond link D_ max=400. We set J as the reference energy scale, and consider J=1 in the following.
As we have done in the CIM, to characterize the two phases we focus on the finite-system AF correlation 
function ( 𝒮^z_[1],l)^2 and on the x-axis string operator 𝒪^x_[1],l
[the corresponding order parameters are defined, in the thermodynamic limit, by Eqs. (<ref>) and (<ref>)].
Since we have ( 𝒮^z_[1],l)^2≥0 and 𝒪^x_[1],l≤0 in the XXZ chain, in the following we will consider the absolute value of the string parameter in order to deal with positive quantities.

Before analyzing in detail the phase transition in the presence of disorder, let us briefly discuss the clean XXZ model.
The situation is summarized in Fig. <ref>, which shows
the behavior of the bulk expectation value of the staggered correlator along z (blue data set),
and of the absolute value of the string correlator along x (green data set).
We see that, at a critical value Λ^(L)_c of the anisotropy term, the SO vanishes
and the staggered order starts to take a finite value: as explained before, this is an indication
of the occurrence of the Haldane-Néel phase transition.
The position of the critical point in the thermodynamic limit, Λ^(∞)_c, can be inferred from the finite-size scaling of the crossing point between the two curves; from our simulations at finite L, we estimate Λ^(L)_c=1.17±0.01, which is in agreement with the value Λ^(∞)_c ≃ 1.186… found in the thermodynamic limit <cit.>.





Based on our knowledge on the clean XXZ model, we now focus on the Haldane-Néel phase transition 
in the disordered scenario. 
Since MPS simulations are computationally more demanding and can only afford systems 
with a comparatively small length, we adopted a procedure  slightly different from the one used
for the CIM, in order to estimate the SO correlator [Eq. (<ref>)] and the AF correlator [Eq. (<ref>)]
from the bulk expectation values of the corresponding finite-size correlators.
The two methods coincide in the thermodynamic limit, but the one described here is more appropriate for the smaller
value of system size and number of disorder realizations which we can obtain with DMRG in the XXZ chain because it enables to
minimize the uncertainty in the averages.

Namely, we are interested in computing the bulk expectation values of a given two-point observable, 
of the form 𝒜̂_k,k+l. 
We recall that both the AF correlator, ( 𝒮^α_[1],l)^2, 
and the SO correlator, 𝒪^α_[1],l, can be seen as expectation values of observables 
which live on a given number of sites in between the k-th and the (k+l)-th site.

After fixing the system size and the disorder realization, we first compute a space average 
over different lengths l of the correlator <cit.>, in order to average out space fluctuations.
We discard the sites that are close 
to the two chain ends, thus disregarding boundary effects 
(see Appendix <ref> for details).
Then, we repeat the simulation by varying the configuration of the disorder, and eventually perform 
a second average of such obtained space averages, over different disorder realizations.

The obtained correlators have an uncertainty  (denoted by error bars in Fig. <ref>) which
is estimated by computing the variance of the space averaged correlators for the different realizations of disorder, as detailed in Appendix <ref>.





We first present our numerical analysis for Λ_j uniformly distributed in the interval [0,Λ_ max], 
for all j. In this case, we expect the system to undergo the Haldane-Néel phase transition 
as Λ_ max is varied across some critical value.

In our simulations we see that, in the Néel phase, the AF pattern is affected by the presence 
of kinks (domain walls, where the AF pattern is reversed), which hide the presence of long-range AF order (<ref>). 
As detailed in Appendix <ref>, the presence of such kinks is a numerical artifact 
due to the non-perfect convergence of the MPS algorithm. Thus, instead of computing the staggered correlator 
as in Eq. (<ref>), we can get rid of the kinks and reveal the presence of AF long-range order 
by computing the bulk average [Eq. (<ref>)] of the Néel correlator, which is defined as


    𝒩^z_[1],l = | ⟨Ŝ^z_k Ŝ^z_k+l⟩|   .


Notice that ( 𝒮^z_[1],l)^2 coincides with the Néel correlator 
in Eq. (<ref>) in the case of perfect AF order (no kinks), 
but differently from the staggered correlator, 𝒩^z_[1],l is insensitive 
to such numerical artifacts, because of the presence of the absolute value in Eq. (<ref>). 
Thus, we characterize Haldane and Néel phases in the disordered chain by looking respectively
at the SO along x [see Eq. (<ref>)] and the Néel correlator 
in the limit of l →∞: 𝒩^z_[1] = lim_l→∞𝒩^z_[1],l.





In Fig. <ref>, we show the result of simulations with L=120 and Λ_ min=0,
after averaging over space and over disorder. 
Blue points correspond to the Néel order, 𝒩^z_[1], 
whereas green points are the SO data, |𝒪^x_[1]|. 
We observe that the Néel order is zero for sufficiently small Λ_ max, and starts to increase 
around a given value of Λ_ max. Conversely, the SO along x is nonzero for small Λ_ max, 
and goes to zero as Λ_ max is increased. This behavior is analogous to the one for the
Haldane-Néel phase transition in the clean XXZ model (Fig. <ref>). 
Furthermore, we can estimate a critical point Λ^(L)_ max,c which is shifted 
with respect to the clean value: for the simulation in Fig. <ref>, 
we find Λ_ max,c^(L)=2.32±0.02. For a given finite size the transition behaves as a crossover, 
exactly as it occurs for the CIM (upper panel of Fig. <ref>) and we estimate the finite-L approximation of the critical
point as the abscissa of the crossing point of the two curves of the correlators. The error bars of the crossing points are estimated as follows. We consider the plot of the correlators (for instance Fig. <ref>) and, being interested in the fluctuations of the intersection of the disorder
averages, we divide the error bars by √(N_ av), according to the central limit theorem. Later, we proceed in a way similar to Fig. <ref>: each time we add or subtract these fluctuations to the disorder-averaged 

correlators and then we  take  the half-dispersion of the 4 resulting values of the crossing.  This gives the uncertainty of the crossing point.


As in the CIM case, in order to extrapolate the value of the critical point 
in the thermodynamic limit, we need to perform a finite-size scaling analysis and
repeat the same simulations as in Fig. <ref> for different values of L. 
The result is shown in Fig. <ref>:
blue points correspond to the estimated values of Λ_ max,c^(L), and uncertainties
are computed as explained above. 
To find the thermodynamic value of the critical point, we show the data as a function of 1/L 
and perform a best fit with the function f(L)=a/L+Λ_ max,c^(∞). 
From the result in Fig. <ref>, 
we extrapolate the critical value in the asymptotic L→∞ limit: Λ_ max,c^(∞)=2.40±0.01.





We can also define the height of the crossing point 
w^(L), as in Sec. <ref>. For each value of L, the uncertainty on the value of w^(L) is computed from that of the disorder-averaged correlators. We show the result in Fig. <ref>, 
from which we see that w^(L) decreases as L is increased. From the fit, we estimate w^(∞)=0.06±0.01. 
This is an indication of the fact that the phase transition becomes sharper and sharper as L is increased.

So far, we have discussed the Haldane-Néel phase transition for Λ_j uniformly distributed 
in the interval [Λ_ min,Λ_ max], for all j, using Λ_ min=0. 
In order to see how the position of the critical point is affected by the choice of Λ_ min 
and Λ_ max, we simulate the disordered model of Eq. (<ref>) 
varying Λ_ max and using Λ_ min=Λ_ max/(n+1), where n is a positive integer number.

Our results for L=120 are shown in Fig. <ref>. 
For each value of n, we estimate the position of the critical point Λ_ max,c^(L)(n) 
as explained for the data in Fig. <ref>. 
If we define the algebraic average of {Λ_j} in the chain, i.e.


    Λ̅=Λ_ min+Λ_ max/2=Λ_ max/2 n+2/n+1  ,


we find that the Haldane-Néel phase transition in the disordered chain occurs when Λ_ max is such that 
the mean value of {Λ_j} in Eq. (<ref>) equals the critical value Λ_c 
of the clean chain, i.e., inverting Eq. (<ref>) and showing explicitly the dependence on L:


    Λ_ max,c^(L)(n)=2Λ^(L)_c n+1/n+2  .


As is evident from Fig. <ref>, the position of the critical point agrees with 
the scaling given by Eq. (<ref>), where we use Λ^(L)_c as fit parameter. 
From the fit, we estimate Λ^(L)_c=1.157±0.002, which is not in disagreement with the clean value Λ^(L)_c ≃ 1.17 
found in Fig. <ref> for L=240. 
We ascribe the slight discrepancy of the two estimates to finite-size effects (we use L=120 for the data in Fig. <ref>).

Differently from the CIM, the study of the thermodynamic singularities in the case of the XXZ model
  requires a much larger computational effort, due to the increased numerical complexity. This is left as an open issue for a future work.









§ CONCLUSIONS AND PERSPECTIVES
 




In conclusion, in this work we have studied the existence of the non-local string order in disordered spin chains.
We have focused on two models, the spin-1/2 cluster Ising chain and the spin-1 XXZ Heisenberg chain, which are
well known to show a transition from antiferromagnetism to string order in the clean case (the first model is moreover interesting for applications
in quantum information <cit.>). We have discovered that this
transition persists in both cases if disorder is added. 

In the disordered cluster Ising model we have found a transition from antiferromagnetism to string order by numerically studying the order parameters in
the ground state; we did this using the Jordan-Wigner mapping on an integrable free-fermion model. We have seen also that the transition manifests
in the fermionic representation as a change in the number of 0-energy edge modes. In this model we have found analytically the position of the
transition point using the strongest coupling renormalization group: this analytical prediction is fully confirmed by the finite-size-scaling
on our numerical results. Moreover, studying the thermodynamic singularity at vanishing temperature, we have found a transition between a non-singular and a singular behaviour when the strength of the disorder is increased. We have seen that this transition point is different from the one separating the SO and the AF phase.

In the disordered spin-1 XXZ Heisenberg chain we have studied the order parameters in the ground state by means of the DMRG technique: we have found a transition between an
antiferromagnetic and a string-ordered phase and we have determined its position by means of finite-size scaling. This model is very interesting
because its string-ordered phase is adiabatically connected to the celebrated Haldane phase <cit.> and it can be experimentally studied thanks
to the new cold-atom techniques <cit.>.

Perspectives of future work include first of all the application of the spin-1 strongest coupling renormalization group <cit.> to the XXZ Heisenberg chain, in order to
analytically predict the phase transition point that we find here numerically. Here
we have only addressed the properties of the ground state: it will be interesting 
to consider the properties of all the spectrum, in connection with many-body 
localization (MBL) <cit.> of interacting non-integrable system. MBL systems can show topological order 
in a large fraction of the excited energy eigenstates <cit.>: it would be interesting to see
if also our non-integrable spin-1 disordered XXZ model shows MBL and if string and antiferromagnetic order persist in excited states. 
A possibility to study these phenomena is applying to the Hamiltonian a quantum quench and look at the dynamics 
of the string correlator. In a clean spin-1 XXZ Heisenberg model the string thermalizes <cit.>, 
only the ground state being ordered, but in disordered systems the situation 
could be much different thanks to the MBL. 







§ ACKNOWLEDGEMENTS

We are grateful to A. Hamma and L. Mazza for fruitful discussions. We thank E. G. Dalla Torre for useful comments
on the manuscript. We acknowledge the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support. M. C. S. acknowledges partial support from the Israel Science Foundation, grants No. 231/14 and 1452/14. A. R. acknowledges financial support from EU through project QUIC, from “Progetti interni - Scuola Normale Superiore” and from his parents. R. F. kindly acknowledges support from EU through project QUIC (Grant Agreement No. 641122), the National
Research Foundation of Singapore (CRP - QSYNC) and the Oxford Martin School.







§ DETAILS ON THE RENORMALIZATION GROUP FOR THE CIM



In the renormalization procedure of the CIM described in Sec. <ref>, we can distinguish
two different types of RG steps: i) the largest coupling is a J_j 
or ii) the largest coupling is a λ_j.

Case i):
If we assume that the largest coupling is J_j, it is not difficult to see that the local Hamiltonian Ĥ_0
of Eq. (<ref>) has four degenerate ground states:


    |1⟩   =   |(+)_j-1 ↑_j (+)_j+1⟩ ,
    |2⟩   =   |(-)_j-1 ↑_j (-)_j+1⟩ ,
    |3⟩   =   |(-)_j-1 ↓_j (+)_j+1⟩ ,
    |4⟩   =   |(+)_j-1 ↓_j (-)_j+1⟩ ,


where |↑_l⟩, |↓_l⟩ are the eigenstates of σ̂_l^z and |(+)_l⟩, |(-)_l⟩ are the eigenstates of σ̂_l^x.

Applying the degenerate perturbation theory to such four ground states, we get the following corrections at first
order in the perturbation V̂ of Eq. (<ref>):



    |ψ_g1⟩    =    |(+)_j-1 ↑_j (+)_j+1⟩
                           +i/2J_j[λ_j-1σ̂_j-2^y|(-)_j-1 ↑_j (+)_j+1⟩
                           +λ_j+1σ̂_j+2^y|(+)_j-1 ↑_j (-)_j+1⟩]  , 
    
      |ψ_g2⟩    =    |(-)_j-1 ↑_j (-)_j+1⟩
                           -i/2J_j[λ_j-1σ̂_j-2^y|(+)_j-1 ↑_j (-)_j+1⟩
                           +λ_j+1σ̂_j+2^y|(-)_j-1 ↑_j (+)_j+1⟩]  , 
    
      |ψ_g3⟩    =    |(-)_j-1 ↓_j (+)_j+1⟩
                           -i/2J_j[λ_j-1σ̂_j-2^y|(+)_j-1 ↑_j (+)_j+1⟩
                           -λ_j+1σ̂_j+2^y|(-)_j-1 ↑_j (-)_j+1⟩]  , 
    
      |ψ_g4⟩    =    |(+)_j-1 ↓_j (-)_j+1⟩
                           -i/2J_j[-λ_j-1σ̂_j-2^y|(-)_j-1 ↑_j (-)_j+1⟩
                           +λ_j+1σ̂_j+2^y|(+)_j-1 ↑_j (+)_j+1⟩]   .


To apply degenerate perturbation theory, we construct and diagonalize the matrix V_ij=⟨i|V̂|ψ_gj⟩
which is given by:


    𝕍=([          -λ_j-1^2+λ_j+1^22J_j λ_j-1λ_j+1J_jσ̂_j-2^yσ̂_j+2^y                             0                             0; λ_j-1λ_j+1J_jσ̂_j-2^yσ̂_j+2^y          -λ_j-1^2+λ_j+1^22J_j                             0                             0;                             0                             0          -λ_j-1^2+λ_j+1^22J_j λ_j-1λ_j+1J_jσ̂_j-2^yσ̂_j+2^y;                             0                             0 λ_j-1λ_j+1J_jσ̂_j-2^yσ̂_j+2^y          -λ_j-1^2+λ_j+1^22J_j ])  .



This matrix can be written as

    V_ij=∑_n^excited states⟨i|V̂|n⟩⟨n|V̂|j⟩/E_ GS-E_n ,

where E_ GS is the energy of the degenerate ground states |i⟩, |j⟩.
Diagonalizing this matrix, one finds the perturbed ground-state eigenenergies at second order in λ: 
 -λ_j-1^2+λ_j+1^22J_j±λ_j-1λ_j+1J_jσ̂_j-2^yσ̂_j+2^y.
We select one of these four eigenstates discarding the others. We arbitrarily choose one of the two states with eigenvalue 
 -λ_j-1^2+λ_j+1^22J_j-λ_j-1λ_j+1J_jσ̂_j-2^yσ̂_j+2^y.
Through the renormalization we have indeed eliminated the site j and generated a new coupling


    -λ_jσ̂_j-2^yσ̂_j+2^y   with  λ_j≃λ_j-1λ_j+1/J_j .


These operators σ̂_j-2/j+2^y are in principle different from the unrenormalized ones: they coincide with them up to terms quartic in λ/J.


Case ii):
If the largest coupling is one of the λ_j, we can reduce to the first case by applying to the Hamiltonian
Eq. (<ref>) the duality transformation <cit.>:


    μ̂_j^x = ∏_k=1^jσ̂_k^z  ,     μ̂_j^z = σ̂_j^xσ̂_j+1^x .


We find the Hamiltonian in the dual representation as


    Ĥ̃̂= - ∑_j [ J_j μ̂_j-1^yμ̂_j^y + λ_jμ̂_j-1^xμ̂_j^zμ̂_j+1^x ] .


Indeed we can see that, in the limit L→∞ we are considering, this Hamiltonian is equal to its dual in Eq. (<ref>),
with λ_j and J_j exchanged. The term with the largest coupling which has to be renormalized
is indeed Ĥ̃̂_0 = -λ_jμ̂_j-1^xμ̂_j^zμ̂_j+^x.
Applying to it the same analysis of the first case, we see that the renormalization procedure eliminates the site j
in the dual representation 
and generates the term in Eq. (<ref>).

It is now easy to show that, after many RG steps, the couplings are renormalized according to


    J_j   =   J_j-2lJ_j-2l+2⋯ J_j+2l-2J_j+2l/λ_j-2l-1λ_j-2l+1⋯λ_j+2l-5λ_j+2l-3 ,
    λ_j   =   λ_j-2lλ_j-2l+2⋯λ_j+2l-2λ_j+2l/J_j-2l+1J_j-2l+3⋯ J_j+2l-3J_j+2l-1 .


Applying the central limit theorem we find


    logJ_j     =     2l(  log J-logλ ) 
           + √(2l)(√( Var[log J] +  Var[logλ])) u_J   ,
    logλ_j     =     2l(logλ-log J) 
           + √(2l)(√( Var[log J] +  Var[logλ])) u_λ  ,


where u_J and u_λ are normally distributed random variables; the averages (…) 
and the variances Var[ … ] are performed over the distributions of J_j and λ_j. 
We thus see that, in the limit of infinite RG steps, Eqs. (<ref>) hold.




§ DETAILS ON THE NUMERICAL ANALYSIS OF DISORDERED HEISENBERG CHAINS



In this appendix, we provide details on the strategy that we adopted in order to compute 
the bulk expectation values of generic two-point observables of the form 𝒜̂_k,k+l,
for the numerical results that have been obtained with the MPS-based algorithm on the spin-1 XXZ Heisenberg chain.
We also comment on the analysis of the presence of domain walls in our simulations. 



 §.§ Bulk expectation values


To compute the bulk expectation values of 𝒜̂_k,k+l,
we first fix L, Λ_ min, Λ_ max, and a given instance of disorder.
For each realization, we numerically compute the space bulk-average by discarding 
a certain number Δ L of sites that are close to the chain ends.
Moreover we consider distances l > Δ L such that, provided we are
far from the transition point, they are larger than the system's correlation length <cit.>.
Near the critical point the correlation length tends to diverge, therefore we have always finite-size effects:
in order to understand the properties of the transition in the thermodynamic limit it is thus very important to perform a finite-size scaling as we do in the main text.

Specifically, if the sites are labeled from 1 to L, 
we choose k=Δ L=0.2 L and average the expectation values from l=l_1=0.3 L 
to l=l_2=0.6 L, i.e.:


    𝒜_ avg,h =
    1/l_2-l_1∑_l=l_1^l_2⟨𝒜̂_Δ L,Δ L+l⟩_h   ,


where the subscript _ avg,h denotes the space average for the h-th realization
of disorder.  Then we repeat the simulation by varying the configuration of the disorder in the chain, 
and perform an average over all the N_ av realizations:


    𝒜_ avg =
      1/N_ av∑_h=1^N_ av𝒜_ avg,h  .






Because of the presence of the random {Λ_j}, the value of ⟨𝒜̂_Δ L,Δ L+l⟩, 
as a function of l, is expected to fluctuate in space. Thus, the expectation value 
in Eq. (<ref>) is affected by an uncertainty, which we estimate via the standard deviation


    σ_𝒜_h^2 = 1/l_2-l_1∑_l=l_1^l_2(⟨𝒜̂_Δ L,Δ L+l⟩_h - 𝒜_ avg,h)^2  .









The fluctuations over the disorder realizations are in turn computed via the variance computed from the {𝒜_ avg,h} in Eq. (<ref>):

    σ^2_𝒜̅ = 1/N_ av∑_h=1^N_ av(𝒜_ avg,h -  𝒜_ avg )^2  .

Finally, to motivate the choice of the spatial averages as in Eq. (<ref>),
we show the behavior of the disorder-averaged string correlator, |𝒪^x_[1],l|,
and the disorder-averaged Néel correlator, 𝒩^z_[1],l,
in Fig. <ref>.
Instead of averaging over space for a given realization of disorder, as discussed before, here we average each value
of |𝒪^x_[1],l|(h) and of 𝒩^z_[1],l(h), for fixed k and l, over N_ av realization
of disorder, where the symbol “(h)” indicates that we are computing the expectation values for the h-th realization
of disorder [see Eqs. (<ref>) and (<ref>)]. Explicitly:



    |𝒪^x_[1],l|     =1/N_ av∑_h=1^N_ av|𝒪^x_[1],l|(h)
    𝒩^z_[1],l    =1/N_ av∑_h=1^N_ av𝒩^z_[1],l(h)    .




For the plots in Fig. <ref>, we choose k=0.2 L, L=120, as in Eq. (<ref>). The red vertical lines limit the interval of l over which we compute the spatial averages in Eq. (<ref>). We stress that, as for the computation of the bulk expectation values, the two ways of averaging (average over space/disorder and then average over disorder/space) are actually equivalent, but performing the disorder-average for each value of k and l, as in Eqs. (<ref>), allows us to visualize the average spatial behavior of the string and Néel correlators.

As we see from Fig. <ref>, apart from spatial fluctuations due to the presence of disorder, the behavior of |𝒪^x_[1],l| and 𝒩^z_[1],l
is in agreement with the expected behavior in the two phases (see Sec. <ref>).
We also see that the choice of l_1=0.3 L and l_2=0.6 L, as in Eq. (<ref>), allows us to capture the average bulk expectation value, at least sufficiently far away from the transition point. When we are close to the Haldane-Néel phase transition (e.g., Fig. <ref>, panels with Λ_ max=2.3), for finite L, our numerical results are affected by finite-size effects, and the correct estimation of the critical point can be then performed only by a finite-size scaling, as explained in Sec. <ref>.













 §.§ Disordered XXZ model with symmetry-breaking magnetic field


As we mentioned in Section <ref>, we can expect the string-ordered phase to be present in the disordered XXZ model as long as the 𝔻_2 symmetry is preserved. We now present some numerical data in order to show an example of destruction of the string order in case a symmetry-breaking term is added to the Hamiltonian. Specifically, we simulate the disordered XXZ model in Eq. (<ref>) with the inclusion of a magnetic field along the x-axis:

    Ĥ_B=Ĥ_ XXZ+B_x∑_jŜ^x_j    .

Since the 𝔻_2 symmetry is broken when B_x≠0 in Eq. (<ref>), we expect SO not to be present in the system <cit.>. The result of a simulation with Λ_ min=0, Λ_ max=1.5 and different values of B_x (in units of J) is shown in Fig. <ref>. In order to highlight the different behaviors at long lengths with respect to the B_x=0 case, we use L=180 for B_x=0,0.2 and L=120 for B_x=0.4,0.6. We compute the disorder average of the string correlator |𝒪^y_[1],l| as in Eq. (<ref>). As we see, the addition of B_x≠0 makes the string correlator decay to zero in the l→∞ limit, and no SO is present in the system, in agreement with the general arguments presented in Refs. PerezGarcia_2008,Pollmann_2010.



 §.§ Domain walls


Let us now consider Λ_ min=0 and explicitly focus on a case with SO (i) 
and a case with AF order (ii).

i) In Fig. <ref>, we show the results of a simulation 
with Λ_ max=1.5 (SO phase) and L=120: 
two-point correlator C^z_[1],l≡⟨Ŝ^z_k Ŝ^z_k+l⟩_h (left panel), 
staggered correlator (𝒮^z_[1],l)^2 (central panel), 
and string correlator |𝒪^x_[1],l| (right panel). 
Here all the correlators have been evaluated over a specific realization of disorder, which we term h. 
To avoid boundary effects, we choose Δ L = 0.2 L, L=120 (thus, we fix k=24). 
For this value of Δ_ max, we see that C^z_[1],l oscillates between positive and negative values, 
for sufficiently small l, and it is damped by an exponential decay. 
This is clearly seen in the behavior of (𝒮^z_[1],l)^2, which exponentially goes to zero 
in the bulk of the chain. On the other hand, the expectation value of the string operator 𝒪^x_[1],l
takes a finite value in the bulk of the chain.

ii) We repeat the same analysis as before, but for Λ_ max=5.0 (AF phase). 
The results are shown in Fig. <ref>. For this value of Λ_ max, 
the two-point correlator C^z_[1],l displays an undamped oscillating pattern, signaling the presence of AF order. 
However, we notice that the pattern reverses at l≃29, i.e., where the data display a kink (domain wall). 
The presence of such a kink suggests that the AF order appears only locally (the system tends to form domains). 
In order to see if the presence of domain walls in the pattern of the two-point correlator C^z_[1],l 
is a physical fact or a numerical artifact, we repeat the simulation M times. We fix the  values of L, 
Λ_ min and Λ_ max and the disorder configuration {Λ_j} and in each repetition we vary
the initial random MPS state |Ψ_ in⟩ at the beginning of the MPS algorithm. 
Our purpose is to verify that different initial random states produce different configurations of kinks 
with different GS energies.

To give an example, we show in Fig. <ref> the result of several simulations 
for M different initial random MPS states: { |Ψ_ in(m)⟩}_m = 1 … M. 
For each value of m, we measure the GS energy,
E_ GS(m) = ⟨Ψ_ in(m)|Ĥ_ XXZ |Ψ_ in(m)⟩,
where Ĥ_ XXZ is the Hamiltonian in Eq. (<ref>), and the number of kinks N_ kinks(m), 
which is obtained from the spatial pattern of ⟨Ŝ^z_j⟩. 
We define the quantity δ E_ GS(m) = [E_ GS(m)- min_m{E_ GS(m)}]/| min_m{E_ GS(m)}|, 
and compare the values of δ E_ GS(m) with the corresponding number of kinks.





As is evident from Fig. <ref>, the configurations with zero kinks are associated 
to the lowest value of the GS energy. Furthermore, we see that the GS energy tends to be larger 
for those configurations having a larger number of kinks. In the present case, we have three configurations 
with zero kinks (all associated to the same GS energy), one configuration with one and two kinks, 
three configurations with three kinks, and two configurations with four kinks. 
Simulations ending up with the same number of kinks may have different GS energy, 
since the configuration of kinks along the chain varies as well. 
We see therefore that the number of kinks and their spatial configuration depend on the choice 
of the initial random MPS state, and that the minimum energy is obtained with zero kinks:
We conclude that the presence of domain walls in the magnetization pattern along the chain 
is a numerical artifact due to the fact that the variational MPS algorithm does not perfectly converge 
to the global minimum of the energy functional.

As a consequence, we expect the true GS to have no kinks: also in the presence of disorder, 
there is long-range AF order in the Néel phase. This justifies our choice of using 
the Néel correlator [Eq. (<ref>)] to estimate the staggered correlator 
( 𝒮^z_[1],l)^2: the Néel correlator does not see the unphysical kinks.


100

Goldenfeld:book
  N. Goldenfeld,  Lectures on Phase Transitions and Critical Phenomena (Avalon Publishing, 1992).

Sachdev:book
 S. Sachdev, Quantum Phase Transitions (Second Edition)" (Cambridge University Press, 2011).
  
Nobel_lecture
  Topological Phase Transitions and Topological Phases of Matter, compiled by the Class for Physics of the Royal Swedish Academy of Sciences (2016), available at
https://www.nobelprize.org/nobel_prizes/physics/laureates
/2016/advanced-physicsprize2016.pdf.

Bernevig:book
  B. A. Bernevig with T. Hughes, Topological Insulators and Topological Superconductors (Princeton Univ. Press, 2013).

Simons:book
  A. Altland and B. Simons, Condensed Matter Field Theory 2nd edition (Cambridge University Press, 2010).

Haldane_1983
  F. D. M. Haldane,  Phys. Lett. A 93, 464 (1983);  Phys. Rev. Lett. 50, 1153 (1983).

Kennedy_1992
  T. Kennedy and H. Tasaki,  Phys. Rev. B 45, 304 (1992).

Sanctuary_2003
  W. Chen, K. Hida, and B. C. Sanctuary,  Phys. Rev. B 67, 104401 (2003).

DegliEsposti_2004
  C. Degli Esposti Boschi and F. Ortolani,  Eur. Phys. J. B 41, 503 (2004).

Ueda_2008
  H. Ueda, H. Nakano, and K. Kusakabe,  Phys. Rev. B 78, 224402 (2008).

DallaTorre_2006
  E. G. Dalla Torre, E. Berg, and E. Altman,  Phys. Rev. Lett. 97, 260401 (2006).

Berg_2008
  E. Berg, E. G. Dalla Torre, T. Giamarchi, and E. Altman,  Phys. Rev. B 77, 245119 (2008).

Rossini_NJ2012
  D. Rossini and R. Fazio,  New J. Phys. 14, 065012 (2012).

Erez_PRB10
  T. Kitagawa, E. Berg, M. Rudner, and E. Demler,  Phys. Rev. B 82, 235114 (2010).

Emanuele_arXiv17
B. Friedman, A. Rajak, A. Russomanno, E. G. Dalla Torre, (2017), arXiv: cond-mat/1708.03400.

Emanuele_arXiv16
  A. Russomanno, B. Friedman, and E. G. Dalla Torre, (2016), arXiv: cond-mat/1611.00659.

Griffiths
  R. B. Griffiths,  Phys. Rev. Lett. 23, 17 (1969).

Mohit
  M. Randeria, J. P. Sethna, and R. G. Palmer,  Phys. Rev. Lett. 54, 1321 (1985).

Bray
  A. J. Bray,  Phys. Rev. Lett. 59, 586 (1987).

Boechat
  A. Saguia, B. Boechat, and M. A. Continentino,  Phys. Rev. Lett. 89, 117202 (2002).

Yona
  Y. Dubi, Y. Meir, and Y. Avishai,  Nature 449, 876 (2007).

MFisher_PRB
  M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher,  Phys. Rev. B 40, 546 (1989).

Basko_Ann06
  D. M. Basko, I. L. Aleiner, and B. Altshuler,  Ann. Phys. 321, 1126 (2006).

Oganesyan_PRA75
  V.Oganesyan and D. A. Huse,  Phys.Rev. B 75, 155111 (2007).

Dasgupta
  C. Dasgupta and S. K. Ma,  Phys. Rev. B 22, 1305 (1980).

Fisher_1995
  D. S. Fisher,  Phys. Rev. B 51, 6411 (1995).

Igloi_review
  F. Iglói, and C. Monthus, Phys. Rep. 412, 277 (2005).

Igloi
  R. Mèlin, Y. C. Lin, P. Lajkó, H. Rieger, and F. Iglói,  Phys. Rev. B 65, 104415 (2002).

Scalapino
  E. H. Kim, G. Fáth, J. Sólyom, and D. J. Scalapino,  Phys. Rev. B 62, 14965 (2000).

Boechat1
  B. Boechat, A. Saguia and M. A. Continentino, Solid State Commun. 8, 411 (1996).

Lajko_PRB
  P. Lajkó, E. Carlon, H. Rieger, and F. Iglói, Phys. Rev. B 72, 094205 (2005).
    
Schollwock2011
  U. Schöllwock, Ann. Phys. 326, 96-192 (2011).
  
Pollmann_2010
  F. Pollmann, A. M. Turner, E. Berg, and M. Oshikawa,  Phys. Rev. B 81, 064439 (2010).
  
PerezGarcia_2008
  D. Pérez-García, M. M. Wolf, M. Sanz, F. Verstraete, and J. I. Cirac,  Phys. Rev. Lett. 100, 167202 (2008).

Son_2011
 W. Son, L. Amico, R. Fazio, A. Hamma, S. Pascazio, and V. Vedral, Europhys. Lett. 95, 50001 (2011).

Smacchia_2011
  P. Smacchia, L. Amico, P. Facchi, R. Fazio, G. Florio, S. Pascazio, and V. Vedral,  Phys. Rev. A 84, 022304 (2011).

Lieb_AP61
  E. Lieb, T. Schultz, and D. Mattis,  Ann. Phys. 16, 407 (1961).

Pfeuty_AP70
  P. Pfeuty,  Ann. Phys. 57, 79 (1970).

Barouch_PRA71
  E. Barouch and B. M. Mc Coy,  Phys. Rev. A 3, 786 (1971).

Pollmann_2012
  F. Pollmann, E. Berg, A. M. Turner, and M. Oshikawa,  Phys. Rev. B 85, 075125 (2012).

DenNijs_1989
  M. den Nijs and K. Rommelse,  Phys. Rev. B 40, 4709 (1989).

Su_2012
  Y. H. Su, S. Y. Cho, B. Li, H.-L. Wang, and H.-Q. Zhou,  J. Phys. Soc. Jpn. 81, 074003 (2012).

Russomanno_JSTAT13
  A. Russomanno, A. Silva, and G. E. Santoro,  J. Stat. Mech. 09012 (2013).

Note_winding
In case of periodic boundary conditions the Hamiltonian of the clean system can be written as the direct sum of 2×2 blocks, each corresponding to a specific value of the momentum k in the one-dimensional Brillouin zone.
The ground state corresponds to a path in the Bloch sphere parameterized by k∈[0,2π]. The winding number corresponds to the Berry
phase <cit.> accumulated in this path, divided by π. Physically, this number counts how many times the Bloch vector turns around a singularity (a point where the gap of the 2× 2 block vanishes), a topological
property which cannot be modified by applying continuous transformations. The only way to modify this property is to make the singularity cross the path,
which implies breaking the continuity of the path (at the crossing point the gap of the 2× 2 Hamiltonian closes and the ground-state Bloch vector is not defined). Therefore, when two phases with different
winding number are put into contact, the gap of the Hamiltonian is forced to close at the boundary, that's why the zero energy boundary modes appear.

Berry
M. V. Berry, Proc. R. Soc. London A 392, 45 (1984).

kitaev2001unpaired
  A. Y. Kitaev,  Phys. Usp. 44, 131 (2001).

Alicea_2012
  J. Alicea,  Rep. Prog. Phys. 75, 076501 (2012).





prb85
Y. Niu, et al., Phys. Rev. B 85, 035110 (2012).

prb88
W. DeGottardi, et al., Phys. Rev. B 88, 165111 (2013).

Note1
  In order to be reliably close to the thermodynamic limit value, one should also ensure that space
  averages are performed over distances l larger than the typical correlation length of the system.
  Therefore, we perform our averages for l> Δ L with Δ L /L = 0.2.
  Near the criticality the correlation length tends to diverge, therefore finite-size effects
  are unavoidable. For this reason, the finite-size scaling we perform is so important
  in order to do statements on the thermodynamic limit behavior.

senko_PRX
 C. Senko, P. Richerme, J. Smith, A. Lee, I. Cohen, A. Retzker and C. Monroe, Phys. Rev. X 5, 021026 (2015).

Chitta_JSTAT13
  B. Bauer and C. Nayak,  J. Stat. Mech. 09005 (2013).
  
Rossini_PRB14
  L. Mazza, D. Rossini, M. Endres, and R. Fazio, Phys. Rev. B 90, 020301(R) (2014).

Strinati_PRB16
  M. Calvanese Strinati, L. Mazza, M. Endres, D. Rossini, and R. Fazio,  Phys. Rev. B 94, 024302 (2016).



