Realising the Symmetry-Protected Haldane Phase in Fermi-Hubbard Ladders

Topology in quantum many-body systems has profoundly changed our understanding of quantum phases of matter. The paradigmatic model that has played an instrumental role in elucidating these effects is the antiferromagnetic spin-1 Haldane chain. Its ground state is a disordered state, with symmetry-protected fourfold-degenerate edge states due to fractional spin excitations. In the bulk, it is characterised by vanishing two-point spin correlations, gapped excitations, and a characteristic non-local order parameter. More recently it was understood that the Haldane chain forms a specific example of a more general classification scheme of symmetry protected topological (SPT) phases of matter that is based on ideas connecting to quantum information and entanglement. Here, we realise such a topological Haldane phase with Fermi-Hubbard ladders in an ultracold-atom quantum simulator. We directly reveal both edge and bulk properties of the system through the use of single-site and particle-resolved measurements as well as non-local correlation functions. Continuously changing the Hubbard interaction strength of the system allows us to investigate the robustness of the phase to charge (density) fluctuations far from the regime of the Heisenberg model employing a novel correlator.

Topology in quantum many-body systems has profoundly changed our understanding of quantum phases of matter. The paradigmatic model that has played an instrumental role in elucidating these effects is the antiferromagnetic spin-1 Haldane chain [1, 2]. Its ground state is a disordered state, with symmetry-protected fourfold-degenerate edge states due to fractional spin excitations. In the bulk, it is characterised by vanishing two-point spin correlations, gapped excitations, and a characteristic non-local order parameter [3,4]. More recently it was understood that the Haldane chain forms a specific example of a more general classification scheme of symmetry protected topological (SPT) phases of matter that is based on ideas connecting to quantum information and entanglement [5][6][7]. Here, we realise such a topological Haldane phase with Fermi-Hubbard ladders in an ultracoldatom quantum simulator. We directly reveal both edge and bulk properties of the system through the use of single-site and particle-resolved measurements as well as non-local correlation functions. Continuously changing the Hubbard interaction strength of the system allows us to investigate the robustness of the phase to charge (density) fluctuations far from the regime of the Heisenberg model employing a novel correlator.
Topological phases of matter often share a deep connection between their bulk and edge properties [8,9]. In the case of the Haldane chain, the bulk exhibits a hidden antiferromagnetic (AFM) order characterised by AFM correlations interlaced with an arbitrary number of S z = 0 elements. This pattern can only be revealed through non-local string correlations that are sensitive to the local spin states, requiring however a detection of the quantum many-body system with microscopic resolution. Even though this was not possible in early experiments on spin-1 chains, evidence for a spin-gap as well as spin-1/2 localised edge states was found using neutron scattering [10,11] or electron resonance experiments [12,13] while not directly probing this hidden order or spatially resolving the edge states. Recent developments in quantum simulations allow one to go beyond such solid-state bulk measurements by observing quantum many-body systems with single-site resolution [14][15][16][17][18] and in a fully spin-and density-resolved way [19,20]. This provides a rich diagnostic tool to obtain a direct microscopic picture of the hidden order in experiments [21,22]. The power of this technique has also been demonstrated recently in a study that was able to reveal a SPT phase in the hardcore boson Su-Schrieffer-Heeger (SSH) model using Rydberg atoms [23]. Here we expand on those results by realising the Haldane phase in a spin system with tunable coupling strength, system size and controlled charge fluctuations. We show this by measuring both topological and trivial string order parameters.
An instructive way to engineer the Haldane phase in systems of spin-1/2 fermions is based on the celebrated AKLT model [4,24], in which a spin-1 particle is formed out of two spin-1/2 particles. Thus, spin-1/2 ladder systems emerge as an experimentally realisable platform for the Haldane phase. While a natural implementation with spin-1 particles on individual rungs requires ferromagnetic rung couplings and antiferromagnetic leg couplings, a much wider variety of couplings in spin-1/2 quantum ladders features the Haldane phase [25,26]. This includes the antiferromagnetic Heisenberg case, which we realise here as the strong-interaction limit of the Fermi-Hubbard model.
In our experiment, we prepare such ladders by adiabatically loading a spin-balanced mixture of the two lowest hyperfine states of 6 Li into an engineered lattice potential (see Methods). As illustrated in Fig. 1a, we realise four isolated two-leg ladders with variable number of L unit cells, surrounded by a low-density bath of particles [27]. The atoms in the lowest band of the optical lattice are well described by the Fermi-Hubbard model with tun-  Figure 1. Probing topological phases in spin-1/2 ladders of cold atoms. a, Realisation of tailored spin-1/2 ladders in a single plane of a 3D optical lattice with a potential shaped by a digital micromirror device (DMD). The dilute wings of the potential are well separated from the homogeneous ladder systems. Using quantum gas microscopy, we obtain fully spin-and density-resolved images of the system. The inset shows a single-shot fluorescence image of the prepared ladder without spin resolution. b and c, Connecting spin-1/2 ladders to integer-spin chains by grouping pairs of spins in unit cells. For diagonal unit cells (b) the AFM Heisenberg ladder adiabatically connects to the Haldane spin-1 chain showing spin-1/2 edge states and hidden long-range order (i.e. AFM order interspersed with S z = 0 unit cells). We realise this topological configuration by blocking one site on each end of the ladder. For the case of straight edges (c), the unit cells coincide with rungs of the ladder and the system is in the topologically trivial configuration. For J J ⊥ , singlets form on the rungs leading to a spin-0 chain. The energy spectra of the systems grouped by total magnetisation M z display gapped fourfold near-degenerate ground states for the topological configuration and a single ground state for the trivial one. Sketch for L = 7.
nelling energies t (chain), t ⊥ (rung) and on-site interactions U . For half-filling and at strong U/t ,⊥ ∼ 13, used throughout most of our experiments (see Methods for details), density fluctuations are suppressed and the spin ladder realises the Heisenberg model [28] with positive leg and rung couplings, J ,⊥ = 4t 2 ,⊥ /U and the spin-1/2 operatorsŜ x,y at site (x, y) with A, B denoting the two legs of the ladder.
The topological properties are most easily explained in the limit J ⊥ J , where strong spin singlets form along the rungs and the system exhibits an energy gap of J ⊥ . The behaviour on the edges of the ladder then depends on how the system is terminated. For tilted edges (see Fig. 1b), two unpaired spin-1/2s remain and the many-body system has a fourfold degeneracy that is only weakly lifted by an edge-to-edge coupling which vanishes exponentially with system size (see SI). In the trivial case of straight edges (see Fig. 1c), all spins pair into singlets and the ground state is unique. These descriptions remain valid even for weaker J ⊥ /J , where the singlet alignment may change between vertical and horizontal, but any line between two rungs cuts an even number of singlets [29,30].
To make the analogy between the spin-1/2 ladder and the Haldane integer chain more apparent, we switch to a description in terms of total spin per kth unit cell, S k =Ŝ k,A +Ŝ k,B , where the indices (A,B) indicate the two spin-1/2s in the same unit cell makingŜ k an integer spin. In the diagonal unit cell such a system shows a high (≥ 80%) triplet fraction [26] (see SI). We note that this spin ladder can be adiabatically connected to a spin-1 chain by including ferromagnetic couplings within the unit cell [25]. However, having a high triplet fraction is not essential for having a well-defined Haldane phase as both systems share the same universal SPT features [26].
The defining property of the Haldane SPT phase is that it is an integer spin chain with spin-1/2 edge modes: the bulk SO(3) symmetry is said to fractionalise into SU (2) symmetry at the edge. It has no spontaneous symmetry breaking and thus the spin correlation function Ŝ z kŜ z k+d is short-ranged. Instead, the aforementioned symmetry fractionalisation [6, 7] can be detected in the bulk using string order parameters [3, 31] with an on-site symmetryÛ l and endpoint operatorÔ k where l denotes the unit cell and d the string distance (see Fig. 2 and SI). This correlator probes the transformation behaviour of the bulk under a symmetryÛ l , e.g. a spin rotation around the z-axis by π,R z l ≡ exp iπŜ z l . The pure-string correlator g 1,R z (d), whereÔ k = 1 and U l =R z l , is non-zero for d 1 if the edge does not have half-integer spins [31]. This is the case for the topologically trivial configuration but not for the Haldane phase The atomic density distribution n of straight-and tiltededge ladders. b, The amplitudes of the spin-string correlator gSz,Rz (green circles) and the string-only correlator g1,Rz (grey squares) observed as a function of the spin distance over d unit cells. The cartoon illustrates the unit cells, the spin total spin S z per unit cell, and the string correlators for a subsystem with d = 3. In the trivial configuration (rung unit cells), |g1,Rz (d)| is well above zero, whereas |gSz,Rz (d)| is rapidly vanishing at d > 1. In contrast, for the topological configuration (diagonal unit cells), |gSz,Rz (d)| shows a long-range correlation, while |g1,Rz (d)| is close to zero. In both cases, the two-point spin-spin correlation C(d) decays rapidly to zero as a function of the distance d (insets). The correlators g1,Rz , gSz,Rz and C(d) are evaluated for fixed total magnetisation M z = 0. c, Amplitudes of the rung-and inversion-averaged local magnetisations |m z (x)| plotted as a function of position x along the chains for different M z . In the imbalanced spin sector of the topological configuration (M z = ±1), the result displays a localisation of the excess spins at the edges signalling the presence of edge states. All data was taken with J ⊥ /J = 1.3(2). Error bars denote one standard error of the mean (s.e.m) and are smaller than their marker size if not visible.
where the symmetry is fractionalised. The spin-string operator g S z ,R z (d) [3],Ô k =Ŝ z k , exhibits the opposite behaviour and is non-zero only in the Haldane phase (see SI for details about the symmetries of the Haldane phase). Thus we can identify the Haldane phase by comparing the two string correlators g S z ,R z and g 1,R z and observe opposite behaviour in the two different regimes.. Another perspective on g S z ,R z can be gained by recognising it as a normal two-point correlator at distance d which ignores all spin-0 contributions along the way ("squeezed space" [22,32]). In the underlying spin-1/2 system, this order stems from N − 1 consecutive rungs dominantly consisting of N −1 singlets and two spin-1/2s which have a combined total spin of +1, 0, or −1.
In order to observe the characteristics of the SPT phase, we prepare a two-leg ladder of length L = 7 and J ⊥ /J = 1.3(2) in both the topological and the trivial configuration. The tailored potential yields a homogeneous filling of the system with sharp boundaries (see Fig. 2a) characterised by a remaining density variance over the system of 2×10 −4 . To focus on the spin physics, we select realisations with N ↑ + N ↓ = 2L per ladder. Additionally, we exclude ladders with an excessive number of doublon-hole fluctuations and do not consider strings with odd atom numbers in the string or the endpoints of the correlator (see Methods). We characterise the spinbalanced ladder systems (M z ≡ (N ↑ − N ↓ )/2 = 0) by evaluation of the string order parameters as defined in Eq. (2). In the topological configuration, we observe a fast decay of g S z ,R z over a distance of ∼1 site and a long-range correlation up to d = 6 with a final value of g S z ,R z 0.1 (see Fig. 2b). In contrast, for the trivial configuration, the correlation decays rapidly to zero as a function of the string correlator length. The opposite behaviour is seen for g 1,R z (d), demonstrating the hidden correlations expected for both phases.
Furthermore, the two-point spin correlation, C(d) ≡ g S z ,1 (d) = Ŝ z kŜ z k+d , yields only the short-range AFM correlation characteristic for a gapped phase (see insets). It is interesting to note that at the largest distance in the topological case, C(d = 6) displays a clear (negative) correlation between the two edge spins, despite small correlations at shorter distances. This (classical) correlation confirms the existence of a non-magnetised bulk, such that spins on the edges of the system must be of opposite direction at global M z = 0.
We probe the edges explicitly by measuring the amplitude of the local rung-averaged magnetisation m z (x) as a function of rung position x for different sectors of the ladder magnetisation M z (see Fig. 2c). In the case of an imbalanced spin mixture with M z = ±1, we see that the two end sites exhibit on average a higher magnetisation than the bulk rungs in the topological configuration. This is consistent with the bulk of the ground states of both phases forming a global singlet and only the edges of the topological phase carrying an excess spin-1/2 without energy cost. The measured bulk magnetisation can be attributed to finite temperature effects (see SI).
The SPT phase is expected to be robust [26] upon variation of the ratio J ⊥ /J , but maintains a finite gap in the system. We realise both the trivial and topological configuration with different t ⊥ /t at almost fixed U and study the string correlators at maximal distance (L − 1) versus J ⊥ /J (see Fig. 3a). For the topological configuration, we observe g 1,R z (L − 1) 0 and |g S z ,R z | > 0 for all J ⊥ /J with a maximum around J ⊥ /J 1.3(2) (i.e. J ⊥ /(J ⊥ + J ) 0.56(4)), while for the trivial case, the role of the correlators is reversed. Both phases continuously connect in the limit of two disconnected chains at J ⊥ = 0. These observations demonstrate qualitatively all the key predictions of the antiferromagnetic spin-1/2 ladder at T = 0 [26] and the strength of the measured correlations are consistent with exact diagonalisation (ED) calculations using an entropy per particle S/N = (0.3 − 0.45) k B (shaded lines in Fig. 3a).
We reveal these features despite the finite temperature in our system, which would destroy the long-range hidden order in an infinite system. The total entropy in our system is, however, still low enough to yield a large fraction of realisations of the topological ground state. In larger systems, the total number of thermal excitations grows (at fixed entropy per particle) and the non-local correlator |g S z ,R z (L − 1)| decreases (see inset of Fig. 3a), consistent with vanishing correlations in the thermodynamic limit thus yielding a restriction on our system size at our level of experimental precision and entropy per particle (see SI). Finite size effects are explored in detail in the SI. We note that the difference between the SPT phase and the trivial phase is here clearly shown by considering both g S z ,R z and g 1,R z .
To investigate the localisation length of the edge states, we evaluate our data for M z = ±1 and plot the local magnetisation per unit cell m z (k) for different system sizes (see Fig. 3b). Due to the singlets in the bulk, the excess spin is most likely found at the edges of the system. This spin partly polarises the neighbouring sites antiferromagnetically leading to an exponentially localised net magnetisation with AFM substructure [33]. The data is well described by the fit function m z (k) = m B + m E (−1) k e −k/ξ + (−1) L−k−1 e −(L−k−1)/ξ with free bulk magnetisation m B , edge magnetisation m E , and decay length ξ. In Fig. 3c, we show how this localisation length ξ decreases as we approach the limit of rung singlets J ⊥ J . Comparison with ED lets us identify two regimes: at J ⊥ J , the measured ξ drops with larger The hidden SPT order is preserved even at low Hubbard interactions as revealed by the novel string correlators g S z ,P ↓ (d) (green circles) and g 1,P ↓ (d) (grey squares) based on the spin-down parity P ↓ . g S z ,P ↓ stays non-zero while g 1,P ↓ is consistent with zero for d = L − 1 over the measured interaction range. The same qualitative behaviour is seen in zero temperature DMRG calculations (shaded line) with L → ∞. c, Spatial distribution of excess magnetisation (M z = ±1) for decreasing U/t . Even far away from the Heisenberg regime, the edge state signal remains strong and only diminishes for very weak U/t . d, Map of zero temperature DMRG (L → ∞) results for the spin string correlator in the entire parameter space of the topological phase. It shows a strictly non-zero g S z ,P ↓ while g 1,P ↓ (L − 1) = 0 everywhere in this phase. The black circles indicate the parameters of the measurements. All experimental data were taken at J ⊥ /J = 1.3(2) and L = 5 in the tilted geometry. M z = 0 in a, b, and d. Error bars denote one standard error of the mean (s.e.m) and are smaller than their marker size if not visible.
J ⊥ and coincides with theory independent of temperature, while at low J thermal effects dominate, limiting the increase of ξ to 3 sites for our system (see SI). Thus far, we have worked in the Mott limit where density fluctuations can be ignored, such that the spin Hamiltonian Eq. (1) is a good effective description of the Fermi-Hubbard ladder. However, it is known that the Haldane SPT phase can be unstable to density fluctuations [34][35][36]. By reducing U/t , the symmetry in the unit cell in the bulk changes from SO(3) to SU (2) as it now may contain both half-integer and integer total spin. This effectively removes the distinction between bulk and edge (see SI). This means that the edge mode and string order parameter are no longer topologically non-trivial, which is also manifested in the fact that the two phases can be adiabatically connected by tuning through a low-U/t regime if one breaks additional symmetries but preserves spin-rotation symmetry [34][35][36]. In particular, the above string orders lose their distinguishing power: g S z ,R z and g 1,R z will both generically have long-range order away from the Mott limit [34].
Intriguingly, despite the breakdown of the above symmetry argument and string order parameter, the Hubbard ladder (with diagonal unit cell) remains a nontrivial SPT phase due to its sublattice symmetry. This symmetry is a direct consequence of the ladder being bipartite (see SI for details). It is simplest to see that this protects the SPT phase in the limit U = 0, where the two spin species decouple, such that our model reduces to two copies of the SSH chain [37]. It is known that such a stack remains in a non-trivial SPT phase in the presence of interactions, i.e. U = 0 [38]. Moreover, together with the parity symmetry of spin-down particles, , it then gives rise to a different string order parameter: the topological phase is characterised by long-range order in g S z ,P ↓ whereas it has vanishing correlations for g 1,P ↓ , with the roles being reversed in the trivial phase. This novel string order is derived in the SI. Remarkably, in the Heisenberg limit, it coincides with the conventional string order parameter used before.
In the topological phase it is meaningful to normalise which effectively excludes endpoints with spin S z = 0. We indeed find unchanged string correlationsg S z ,P ↓ and g 1,P ↓ down to the lowest experimentally explored value U/t = 2.5(2) (see Fig. 4a,b) and edge state signals down to U/t = 5.0(3) (see Fig. 4c). DMRG calculations for L → ∞ confirm non-zerog S z ,P ↓ (L − 1) at T = 0 and for all rung coupling strengths (see Fig. 4d), while g 1,P ↓ (L − 1) is strictly zero. Due to the normalisationg S z ,P ↓ (L−1) goes to 1 for J ⊥ J . In our work, we realised a finite temperature version of the topological Haldane SPT phase using the full spin and density resolution of our Fermi quantum gas microscope. We demonstrated the robustness of the edge states and the hidden order of this SPT phase in both the Heisenberg and the Fermi-Hubbard regime. In the future, studies may extend the two-leg ladder to a varying number of legs where one expects clear differences between even and odd numbers of legs [39] and topological effects away from half-filling [40], or investigate topological phases in higher dimensions [41]. Furthermore, the ladder geometry holds the potential to reveal hole-hole pairing [42] at temperatures more favourable than in a full 2d system. Competing interests: The authors declare no competing interests.
Data availability: The datasets generated and analysed during the current study are available from the corresponding author on reasonable request.  In each experimental run, we prepare a cold atomic cloud of 6 Li in a balanced mixture of the lowest two hyperfine states (F = 1/2, m F = ±1/2). For evaporation, we confine the cloud in a single layer of a staggered optical superlattice along the z-direction with spacings a s = 3 µm and a l = 6 µm and depths V s = 45 E s R and V l = 110 E l R , where E R denotes the recoil energy of the respective lattice. The atoms are harmonically confined in the xy−plane and the evaporation is performed by ramping up a magnetic gradient along the y-direction (see [20]). The final atom number is tuned via the evaporation parameters.
The cloud is then loaded into an optical lattice in the xy−plane with spacings a x = 1.18 µm and a y = 1.15 µm, which is ramped up within 100 ms to its final value ranging from 5 E R to 15 E R depending on the chosen Hubbard parameters. The scattering length is tuned from 230 a B during evaporation, with a B being the Bohr radius, to its final value ranging between 241 a B and 1200 a B , using the broad Feshbach resonance of 6 Li. An overview of the parameters of each dataset is given in the Table I. Simultaneously with the lattice loading, a repulsive potential is ramped up which compensates for the harmonic confinement generated by the curvature of the Gaussian lattice beams and divides the resulting flat area into four disconnected ladder systems along the y−direction (see below "potential shaping"). We achieve temperatures of k B T ∼ 0.9(3) J for the parameters in Fig. 2. For detection, the configuration is frozen by ramping the xy−lattices to 43 E xy R within 250 µs. A Stern-Gerlach sequence separates the two spin species into two neighbouring planes of the vertical superlattice, which are then separated to a distance of 21 µm using the charge pumping technique described in [20]. Finally, simultaneous fluorescence images of the two planes are taken using Raman sideband cooling in our dedicated pinning lattice with an imaging time of 2.5 s [43]. The fluorescence of both planes is collected through the same high-resolution objective. The light is then split into two paths using a polarizing beam splitter. One of the beams passes through a variable 1:1 telescope before both paths are recombined on a second polarizing beamsplitter with a small spatial offset. This allows us to image both planes in a single exposure, each plane in focus on a separate fixed position of our camera. We calibrated the relative position on the camera of the two imaged planes using 300 shots of a spin-split Mott insulator and the matching algorithm described in the supplement of [20]. The overall detection fidelity per atom is 96(1) %.  Table I. Experimental parameters. The parameters system size L, leg coupling t , rung coupling t ⊥ , interaction U and the resulting ratio J ⊥ /J are shown for all datasets. The uncertainties are given for J ⊥ /J and originate from a 5% uncertainty on the hopping parameters t ⊥ and t . For the length scan we keep all other parameters constant, whereas the J ⊥ /J scan demands a tuning of both tunnelling amplitudes in order to keep both U/t and U/t ⊥ high. Where the topologically trivial geometry is realized, it has the same parameters as the topological geometry.

Potential shaping
The ladder systems are created by superimposing the optical lattice with a repulsive potential, which is shaped by projecting incoherent light at 650 nm (generated from a SLED by Exalos EXS210030-03) from a digital micromirror device (DMD) through the high-resolution objective. Four ladders are created by blocking lattice sites with a potential V = 3.5(5), U around each ladder. The area outside of the walled ladders is lifted above the inner ladder potential but remains below the interaction energy U . It thus serves as a reservoir for surplus atoms, which occupy this region once the lowest Hubbard band of the ladders is filled. The flatness of the potential is adjusted for each parameter setting since the intensity of the lattice beams directly influences the curvature of the potential. This is accomplished by realizing a system with about 20% doping and returning the average density of 100-150 experimental runs as feedback to the DMD pattern. We repeat the feedback until we reach a sufficiently flat density profile with a variance < 1 × 10 −3 over the 8L ladder sites. To adjust for drifts in the lattice phase, we continuously track the lattice phase of each experimental run and feedback to the potential position accordingly. In Fig. S1, the average density and occupation histograms of all four ladders and the reservoir area are shown for the dataset of L = 7.

Data selection
In each experimental run, four ladder systems are realised. To fulfil the criteria of the Heisenberg regime, we then select on ladder instances with atom number N = 2L and restrict the total magnetisation to M z = 0, |M z | = 1, or |M z | ≤ 1, depending on the observable, and specify the magnetisation sector whenever data points are presented. |M z | ≤ 1 for 87.5 % of all data. Additionally, for all measurements in the Heisenberg regime, we remove ladders with more than two doublons as those indicate a mismatch of the DMD pattern relative to the lattice phase. To give a specific example, we here give the precise numbers for the data presented in Fig. 2. This dataset consists of 7533 realisations with four ladders each. Out of those 28128 ladders, 6721 have an atom number of 14. 2636 ladders then have a total magnetisation M z = 0, 3094 have a magnetisation fo M z = ±1. Finally of those 2636, 77 have more than 2 doublon-hole pairs, which we exclude as those are most likely caused by drifts of the potential pattern given by the DMD. This leaves then in total 2559 ladders out of the initial 28128 for the calculation of the string correlator. For calculating the string correlators g S z ,R z and g 1,R z at fixed d, we exclude realisations with odd atom number in the bulk area (grey area in the cartoon of Fig. 2b) as those would lead to imaginary contributions to the correlators and also exclude odd atom number at the edge areas (green in the cartoon of Fig. 2b). These cases are mostly due to the finite U/t which still allows for some particle fluctuations. We keep other particle-hole fluctuations like those occurring within the string. Those do not alter the observed string correlation relative to the Heisenberg model.

Nearest-neighbour spin correlations
The two-leg ladder systems show strong antiferromagnetic spin correlations whose dominant orientation depends on the ratio of couplings J ⊥ /J and whose strength is measured by C x (d) = 4 Ŝ z i,jŜ the rung coupling J ⊥ , the nearest-neighbour spin correlator C y along the rung almost vanishes, whereas correlations reach C x (1) = −0.500(6) along the leg direction. For a dominating rung coupling J ⊥ , the C y reaches −0.58(1), indicating singlet formation between the two sites of a rung. Fig. S2a shows the measured spin correlations along both rung and leg for different values of J ⊥ /J . The values match the finite temperature Heisenberg model for an entropy of S/N = (0.3 − 0.4) k B per particle obtained from ED.

Theory simulation
In this work, we have employed two different numerical methods to obtain theoretical predictions for the experimentally measured observables. The results in the Heisenberg regime were obtained using exact diagonalisation (ED) of our spin-1/2 ladders up to sizes of L = 9 (limited by computational resources). For each data point, the system size and geometry in the ED simulation is the same as in the experimental data. The finite temperature results were obtained by using the full spectrum. We specify the entropy per particle s = S/N , which we find to be approximately independent of coupling parameters in the experimental realisations. The results in the Hubbard regime are calculated using the Density Matrix Renormalisation Group (DMRG) Ansatz [44] based on the TeNPy library (version 0.3.0) [45]. For all calculations, we conserved the total particle number and the total magnetisation. For the phase diagram in Fig. 4d we used the iDMRG technique to obtain the ground state and the values of the string order parameters in the thermodynamic limit. For this, we evaluated the ground state for each parameter and used a maximal MPS bond dimension χ = 1200. The bond dimension is increased in steps ∆χ = 40 and the simulation stopped when the difference in the ground state energy per unit cell

Triplet fraction in the unit-cell
The mapping of the trivial configuration to a spin-0 chain and of the topological configuration to a spin-1 chain (see Fig. 1 of the main text) is investigated numerically in ED calculations on a system of length L = 5 at zero temperature with M z = 0 (see Fig. S3). As expected, the rung singlet fraction increases monotonically with J ⊥ /J and approaches 1 for J ⊥ J . Remarkably the triplet fraction along the diagonals is always high (≥ 80%) and reaches its maximum at J ⊥ /J ∼ 1 consistent with [26]. The order parameters in the main text prove that, even when tuning J ⊥ /J away from the limiting cases, the system stays within the respective trivial (topological) phase.

Normalisation effects
The string-only correlator g 1,R z is naturally normalised as it returns 1 for any state with an even number of |S z | = 1 in the string. The spin-string correlator g S z ,R z only equals 1 in a classical spin-1 Néel state. Any spin-1 state with rotational symmetry has g S z ,R z < 1 due to the presence of some S z = 0 at the end of the string. A meaningful normalisation is given byg S z ,R z = η g S z ,R z with η −1 = Ŝ z k Ŝ z k+d describing the probability that neither endpoint of the string has spin S z k = 0. In direct analogy to Bayesian conditional probabilities,g S z ,R z describes the string-correlation between |S z | = 1 spins and thus |g S z ,R z | = 1 for the AKLT state.
The application of this normalisation to the data of Fig. 2 is shown in Fig. S4a with the explicit values for η given in Fig. S4b. For d 1, it is equivalent to the normalisation used in [26]. We use the same normalisation in the Hubbard regime (see Fig. 4 of the main text).

Finite system size and temperature
In the thermodynamic limit, the SPT phase only exists at strictly zero temperature (i.e. the spin-string correlator g S z ,R z (d → ∞) = 0). The reason is that there are only four ground states, but infinitely many excited states just above the energy gap resulting in infinitely many singlets that can be broken. In our experimental setup, we can nevertheless observe the characteristics of this phase. The finite length of our system limits the number of low energy states available, such that even at a temperature around the gap energy, the ground state is still largely populated (see Fig. S5). The colder the system, the longer the length at which the ground state still dominates. Fig. S5a shows the ground state population for short system lengths at J ⊥ /J = 1.2 for S/N = 0.3 k B and S/N = 0.45 k B corresponding to a temperature of T = 0.6 J , 1.2 J . The ground state population quickly drops as the number of available states increases. The effect of the reduced ground state occupation can be seen in the measured string correlator g S z ,R z (d) in systems of different length (see Fig. S5b). The system of L = 11 shows a much lower value for the string correlator even for short distances, where smaller systems show significantly higher correlations. This restricts the system size up to which signatures of the Haldane phase can be detected in an experiment but even below this system size, all qualitative features of the zero-temperature phase are already present. Furthermore, to show that the signal is not dominated by the small system size, we use DMRG calculations at infinite length at T = 0.9J J ⊥ /J = 1.3 (shaded line in Fig. S5b).

Finite size offset
In a finite-size system, the string correlator g S z ,R z does not approach zero even when the temperature increases to infinity if the magnetisation is fixed. The lack of free fluctuations of M z introduces correlations in the system even for a random distribution of the spins [22]. These correlations, which can be derived from combinatorial considerations, do not depend on the edge termination or coupling parameters. Fig. S6 investigates the effect of this offset on our measurements. As can be seen in Fig. S6a,b, the offset mostly takes sizeable values at the shortest and longest distances of the system. But even at these points, our signal clearly exceeds the offset. Fig. S6c shows the string correlator as a function of the coupling strength (as in Fig. 3a). Interestingly the string correlation value can coincide with or even be lower than the infinite temperature offset for couplings far away from the symmetric point of J ⊥ ≈ J . This is however not an artefact of the measurement but is also reflected in the zero temperature ED calculations (shaded line). Thus it is not meaningful to simply subtract this offset. Similarly we show the offset for the Hubbard regime in Fig. S6d. As noted before, the offset is considerably smaller for intermediate string lengths. Therefore we investigate the value of the string correlator for d = L − 2 (see Fig. S6e). In this case, the offset is considerably smaller than the measured values for any coupling on the topological side. Therefore our system is not dominated by its finite system size but exhibits the topological properties of the bulk characteristic to the Haldane phase.  Fig. 4). e, String correlator gSz,Rz at distance d = L − 2. For all coupling parameters J ⊥ /J in the topological regime, the measured value clearly exceeds the infinite temperature finite size offset (dashed line).

Unit cell and edge effects
In the main text we present two different edge terminations as realisations of the topological and trivial phase. The influence of the edges does however not fully determine the behaviour of the two string correlators. This is to be expected as the string correlators probe the bulk of the system which one could also calculate in an infinite system. Only in conjunction with a chosen unit cell does it make sense to characterise the system via a specific phase.
As an illustrative example we use a singlet chain (see Fig. S7), where every spin-1/2 particle is paired into a singlet with a fixed neighbour. In the case of periodic boundary conditions (a, b, e), the topology of the system with a two-site unit cell is only set by the choice of unit cell: Unit cells around the the singlets result in a trivial bulk, while the shifted unit cell leads to a topological Haldane phase. We stress that on the level of spin-1/2 particles both states are identical and only the different choices of pairing in the analysis lead to the different topologies. By cutting, the ring can be turned into a chain with edges. But because there is no long-range entanglement in the system, the bulk properties stay unaffected by the cut i.e. the topology is still set by the choice of unit cell and not by the position of the cut. Only one of the two choices of units cell, however, agrees with the cut such that no unit cell is split (c, g). By providing such a natural choice of unit cell, the edge is linked to the topology of the bulk. When only considering the bulk it is also possible to choose the opposite unit cell, which disagrees with the cut, leading to the opposite string-correlation results (d, f).
In the main text, we presented results for the natural choice of units cells (vertical for straight edges, diagonal for tilted edges). In Figure. S7h,i, we compare the string correlators g S z ,R z and g 1,R z with both edge terminations at fixed unit cells. In tilted unit cells (h), we always find a topological bulk, while straight unit cells (i) show the trivial correlations. The physical presence or absence of an unpaired spin-1/2 does not affect the bulk significantly. This is consistent with the discussion in the last paragraph and demonstrates, that, even in our relatively small systems, the properties of the bulk are independent of the edge.

Microscopy of edge states
Due to our microscopic resolution, we can study the magnetisation pattern within the localised edge state in detail. In Fig. S8, the spatial magnetisation distribution of a system with strong rung coupling is compared to the balanced situation (J ⊥ ≈ J ). As the M z = 1 sector has positively polarised edge states at T = 0, the magnetisation maps directly reveal the structure of the state: the excess magnetisation dominantly sits at the edge but leaks into the bulk where it induces a staggered magnetisation pattern close to the edge due to the AFM spin coupling. How far the correlation extends into the bulk is set by the leg coupling J relative to the bulk gap. At finite temperature, there is, in addition, some homogeneous magnetisation of the bulk. In Fig. 3 of the main text, we show the unit-cell average of this data.

Edge state splitting
For finite system length, the four ground states of the SPT phase are not truly degenerate but exhibit a finite energy splitting δ. This energy splitting arises from an exponentially suppressed but non-zero overlap between the edge states, coupling them into singlet and triplet states. Whether the singlet state or the triplet state is lower in energy depends on the parity of the system length L, whereas the energy splitting depends on the system length directly (see Fig. S9a).
Experimentally we cannot observe the splitting directly due to our finite temperature but we explore the underlying spin correlations responsible for the energy splitting: we investigate the effect of even and odd system length by comparing a system of L = 6 (see Fig. S9b, c) to one with L = 5 (see Fig. S8). The staggered magnetisation pattern seen for L = 5 is not visible for L = 6 because the induced spin pattern of both edges are incommensurate with each other leading to the higher energy of the triplet state. The opposite is true in the M z = 0 sector, where we find stronger alternating patterns in the even ladder length compared to the odd one (see Fig. S9c around d = 3, 4). Here we analyse the spin-correlations C(1, d) = 4 Ŝ z i,jŜ z i+1,j+d because for M z = 0 the local magnetisation is zero everywhere. These observations illustrate how the alternation between singlet and triplet ground states with L is linked to the AFM polarisation due to the edges.

Localisation length
We next relate the experimentally extracted edge state decay length ξ to the localisation of the edge at zero temperature. Numerically, the length over which the edge modes delocalise can be readily extracted from the aforementioned energy splitting δ The experimental edge decay length is determined, however, using the staggered magnetisation that arises for where k denotes the position of the unit cell along the chain.
These two approaches are numerically compared in  Figure S9. Edge states for even and odd system length. a, Ground-state energies of the Haldane phase for different system lengths. The lowest-lying state alternates between the singlet and triplet state, depending on the parity of the system length. The data is produced via exact diagonalisation of the spin-1/2-Heisenberg ladder. b, The magnetisation map of our system at total magnetisation |M z | = 1 at L = 6 and J ⊥ /J = 1.3(2). It shows the onset of an alternating magnetisation pattern close to the edge sites, but gets lost in the centre of the ladder since the patterns of the two edges have opposite phase. The plotted value m z sgn is the same as in Fig. S8. c, Spin correlations C(1, d) for a system at total magnetisation M z = 0 of L = 6 (brown) and of L = 5 (grey). Both show a strong nearest-neighbour correlation as well as a strong edge to edge correlation, indicating the edge states of opposite spin for this magnetisation. However, for L = 6 the correlations alternate around a finite-size finitetemperature offset, whereas for L = 5 they do not alternate because an alternation does not match the length and the negative endpoints. and a bond dimension χ = 1000 keeping the maximal energy truncation error below 10 −7 . For the edge magnetisation, we calculated the ground state in the sector M z = 1 for L = 100 and a bond dimension χ = 1000 again with an error below 10 −7 for all parameters. Both quantities agree with deviations of less than 6% of their values, confirming the validity of our method to extract the localisation length from the experimental data (see Fig. 3). We furthermore compute the bulk correlation length, which is a direct measure of the bulk gap. It can easily be obtained from the ground state of the infinite system [46] and coincides with the former length scales for most parameters.
However, the decay length obtained from the staggered magnetisation shows a strong temperature dependence (see Fig. S10b), whereas the edge state splitting is a property of the spectrum, and thus independent of temperature. The staggered magnetisation arises from the antiferromagnetic correlation of the spin, which decreases with temperature. We find that the edge state cannot delocalise beyond the thermal coherence length of the system. The finite temperature decay length thus follows the zero temperature decay length when it is small (large J ⊥ /J ) but saturates in the low J ⊥ /J regime at an upper bound given by the temperature of the system.  Figure S10. Decay length comparison. a, The decay length of several quantities is determined using DMRG for a system of L > 50, and U/t = 13. The purple squares show the decay length of the magnetisation pattern, as it is also evaluated on the experimental data. The green circles show the edge localisation length derived from the edge state splitting and the blue triangles show the bulk correlation length. When markers are not visible they coincide with the purple squares. b, Finite temperature results for a system with L = 5 and S/N = (0 − 0.5) kB (from dark to light purple). The decay length is calculated from the magnetisation pattern in the Heisenberg model using ED.

Symmetry fractionalisation
The key to understanding one-dimensional SPT phases is the notion of symmetry fractionalisation [6, 7]. IfÛ = nÛ n is an on-site unitary symmetry, it can be argued that if the ground state |ψ is symmetric underÛ and has a finite correlation length, then n k=mÛ k |ψ =Û LÛR |ψ , whereÛ L (Û R ) is a unitary operator which is exponentially localised to the left (right) of the block of sites k = m, m + 1, · · · , n − 1, n. Moreover, ifV is another symmetry, then the group relations between the fractionalised symmetriesV L andÛ L are the same as those between the bulk symmetriesV andÛ up to potential phase factors. As an example, suppose thatÛ andV commute and suppose thatÛ L,R are bosonic operators (which is the case for all symmetries considered in this work), then Since the fractionalised operators on the left and right have disjoint support (up to exponentially small overlaps), we have thatÛ LVLÛ −1 LV −1 L = e iα 1, i.e.,Û LVL = e iα V LÛL : the fractionalised symmetries commute up to a phase. More generally, the group relations of the fractionalised symmetries define a projective representation of the symmetry group. Some of these phase factors can be gauged away by redefiningÛ L ,Û R → e iβÛ L , e −iβÛ R , while other phase factors are invariant. The collection of such invariant phase factors define the so-called second group cohomology class H 2 (G, U (1)): any non-zero element in this class represents a non-trivial SPT phase. The Haldane SPT case corresponds to where the bulk symmetry group is SO(3) (or its Z 2 × Z 2 subgroup) and the fractionalised symmetries form SU (2) (or its quaternion subgroup).

From symmetry fractionalisation to edge modes
Note that a non-trivial projective representation for the fractionalised symmetries automatically implies edge modes: for open boundaries, one can consider Eq. (3) as acting on the whole system, which moreover implies that U L andÛ R are genuine symmetries of the ground state. Since a projective representation can never act on a onedimensional Hilbert space, there must be a degenerate zero-energy Hilbert space associated to the boundaries (in other words, the ground state cannot be a simultaneous eigenstate of all fractionalised symmetries).

From symmetry fractionalisation to string order parameters
For an on-site symmetryÛ = nÛ n , consider a string operatorÔ † mÛm+1 · · ·Û n−1Ôn with endpoint operator O n and n − m being much larger than the correlation length. Using symmetry fractionalisation, its expectation value can be expressed as Hence, long-range order (LRO) for this string operator is equivalent to Ô †Û L = 0. In the assumption that the ground state is symmetric (i.e., no spontaneous symmetry breaking), a local operator can only have a non-zero expectation value if it is neutral under the symmetry group. This means that LRO is only possible ifÔ is chosen to have the same symmetry charges as the fractionalised symmetryÛ L (i.e., finding the rightÔ allows to infer the projective representation of the fractionalised symmetries). This is how string order parameters can be used to diagnose an SPT phase [31].
As an example, consider an on-site Z 2 × Z 2 symmetry, such as the π-rotationsR x andR z for an integer spin chain. The fractionalised symmetries either commute or anticommute:R L xR L z = ±R L zR L x ; the anticommuting case gives rise to the topological Haldane phase. In particular, this means thatR xR Hence, a string forR z can only have LRO if the endpoint operator is odd underR x . This explains the well-known string order parameter for the Haldane SPT phase: · · ·R zRzRzŜ z . Conversely, the trivial phase can only have LRO for theR z string if the endpoint is even underR x ; choosing the endpoint operator to be the identity operator does the trick.
If one tunes away from the Mott limit, the distinction between these two cases breaks down. The integer spin chain now does not have SO(3) symmetry but rather SU (2) symmetry. More concretely,R x andR z are no longer Z 2 symmetries: they square to fermion parity sym-metryP , not to the identity. (In the Mott limit, fermion parity is a classical number in each unit cell, giving rise to an effective spin chain.) This means that the fractionalised symmetries now obeyR L xR L z =P LRL zR L x , which allows to adiabatically connect the spin chains wherê P L = ±1 as has been demonstrated before [34][35][36].
A novel string order parameter for anti-unitary symmetry LetT be an anti-unitary symmetry. Analogous to Eq. (3), its symmetry fractionalisation can be written aŝ T =Û LÛRK , whereK is complex conjugation defined with respect to a basis that factorises between left and right (see [6,38] for details). A single anti-unitary Z T 2 symmetry can protect a non-trivial SPT phase. More precisely,T 2 implies thatÛ LKÛLK =Û RKÛRK = e iθ = ±1; the case θ = π is the topological Haldane phase, where the edge mode is a Kramers pair underT .
Usually, it is said that there is no simple string order parameter to detect such an SPT phase protected by Z T 2 . The simplest "string order" is rather involved, requiring two copies of the system and a partial swap [31]. Here, we show that a conventional string order parameter can be constructed if the system has an additional unitary Z 2 symmetry, which we denote asP .
Suppose that we consider phases which are not protected by the combined symmetryPT . We now show that the string order parameter for aP -string, i.e., O † mPm+1 · · ·P n−1Ôn , can be used to diagnose whether the phase is in a trivial or topological phase with respect to the anti-unitary symmetryT . To see this, first note that the phase being trivial with respect toPT implies that ifP =P LPR is the fractionalisation ofP , then PT =P LÛLPRÛRK must obeyP LÛLKPLÛLK = +1. Moreover, sinceP 2 = 1, we can chooseP 2 L = 1, such that We thus obtainTP LT =Û LÛRKPLÛLÛRK = Û LKPLÛLK Û RKÛRK = e iθP L . Hence, whether P L commutes or anticommutes withT encodes what phase we are in. Using the reasoning of the previous section, one can conclude that the string order parameter forP with endpoint operatorÔ (where we chooseÔ to be hermitian) has long-range order ifTÔ nT = e iθÔ n .

Application to Hubbard chain
We have already explained why away from the Mott limit, we can no longer rely on the conventional string order parameter to characterise an SPT phase. However, as pointed out in [38], the (bond-alternating) Hubbard chain is still in a non-trivial SPT phase protected by an anti-unitary Z T 2 symmetry. This uses the bipartite structure of our model (which is evident in the fact that the system lives on a ladder), leading to the symmetry as defined byT : c x,y,s ↔ (−1) x+y c † x,y,s where x = [0, L], y = 0, 1 (corresponding to A or B) and s =↑, ↓. Here the bipartite property is encoded in the factor (−1) x+y as the ladder Hamiltonian only couples sites with opposite parity. This follows directly from combining two facts: (i) it is a Z T 2 symmetry of the model, even away from the Mott limit, and (ii) in the Mott limit, it can be argued that it coincides with spinful time-reversal symmetry [38], which is known to protect the Haldane SPT phase [6].
To construct a string order parameter for this phase, we use the result obtained in the previous section. In particular, consider the additional Z 2 symmetryP ↓ , which is the fermion parity of the down-spin species. In the Mott limit, this symmetry becomes indistinguishable fromR z . From this, we learn thatP ↓T does not protect the SPT phase (which was a condition that we assumed in the previous section). As derived above, this implies that the string operator associated toP ↓ can be used to read off the topological invariant: we are in the Haldane SPT (trivial) phase if the string has long-range order for an endpoint operator that is odd (even) underT . For instance, for the topological phase, we can thus choose the endpoint operatorÔ n =Ŝ z n = 1 2 ĉ † n,↑ĉn,↑ −ĉ † n,↓ĉn,↓ , which is odd under the above anti-unitaryT symmetry.