Operational Entanglement of Symmetry-Protected Topological Edge States

We use an entanglement measure that respects the superselection of particle number to study the non-local properties of symmetry-protected topological edge states. Considering $M$-leg Su-Schrieffer-Heeger (SSH) ladders as an example, we show that the topological properties and the operational entanglement extractable from the boundaries are intimately connected. Topological phases with at least two filled edge states have the potential to realize genuine, non-bipartite, many-body entanglement which can be transferred to a quantum register. We show, furthermore, that the onset of entanglement between the edges can be inferred from local particle number spectroscopy alone and present an experimental protocol to study the breaking of Bell's inequality.

Entanglement in a quantum system is an essential resource for quantum algorithms and for quantum encryption keys. In order to use this resource, an important question is how much of the entanglement present in a quantum state is accessible and can be transfered by local operations to a quantum register consisting of distinguishable qubits. A fundamental difficulty in addressing this question for many-body states of itinerant particles is that the entangled entities itself are indistinguishable. Suppose observers, Alice and Bob, have access to two spatially-separated parts of a quantum system of indistinguishable particles with a conserved particle number N. Then there is always a local operation that will collapse the pure or mixed quantum state ρ AB they share into a state ρ n A ,n B AB with fixed local particle numbers n A (n B ) for Alice (Bob). The operational entanglement (also called accessible entanglement or entanglement of particles) which can be transfered to a quantum register is thus given by [1,2]  (1) where p(n A , n B ) = tr ρ n A ,n B AB is the probability to project onto a state with n A , n B particles in the two subsystems and E[ρ] is an entanglement measure applied to the normalized projected state. The superselection of particle number also gives rise to an additional non-local resource associated with the particle number fluctuations [3,4]. To characterize this second, complementary resource we introduce the generalized number entropy (Shannon entropy) E n = − n A ,n B p(n A , n B ) ln p(n A , n B ) . ( We note that E n = 0 if n A and n B are fixed, E n has an upper bound E max the operational or configurational entropy is now given by Eq. (1) with E[ρ] ≡ E vN [ρ] [5][6][7][8] and the restriction n B = N − n A is placed on the sums in Eqs. (1,2). The number entropy for a bipartition has recently been measured in a cold atomic gas experiment [9] and can be used to obtain a bound on the total entanglement entropy [10,11].
Long-range entanglement in many-body systems is, in general, very fragile against local perturbations. An exception are edge states in systems with an excitation gap in the bulk which exist due to a global symmetry. In this case, the number of edge states is topologically protected and cannot be changed without closing the gap or breaking the symmetry [12,13]. Since topological edge states can exist simultaneously on multiple edges, a filled edge state is a non-local quantum resource.
In this paper we will show that symmetry-protected topological edge states in lattice models can be a source of genuine, spatially-separated, non-bipartite many-body entanglement which can be transfered to a quantum register. We will show, furthermore, that the onset of entanglement between the edges can be inferred from local particle number spectroscopy, making it measurable, for example, by current-day cold atomic gas experiments [9,14,15]. We also present an experimental protocol to test Bell's inequalities [16][17][18][19] and find that for two or more filled edge states sufficient operational entanglement can be extracted from the mixed state of a manybody system of arbitrary length to break them.
In the following, we will consider lattice models of itinerant spinless fermions. Let us first briefly discuss why in this case a single edge state is insufficient to produce any operational entanglement. Let us assume that the edge state is very sharply localized at the left and right edges of a chain which are controlled by Alice and Bob, respectively. To simplify the argument, we ignore the bulk completely for now and assume that the edge state is a single-particle pure state. A maximally entangled edge state in occupation number representation is then of the form |Ψ = (|10 ±|01 )/ √ 2 and has von-Neumann entanglement entropy E vN = ln 2. If Alice and Bob measure their local particle numbers in order to perform local operations, this state however collapses onto the product state |Ψ 1,0 ∝ |10 or |Ψ 0,1 ∝ |01 . Thus, this state has only number entropy E n but no operational entanglement. We therefore need at least two filled edge states to have any operational entanglement. In the latter case, the projected state |Ψ 1,1 can have operational entanglement while the states |Ψ 2,0 and |Ψ 0,2 are again product states.
Model.-In order to build a system with more than one edge state, we couple half-filled Su-Schrieffer-Heeger (SSH) chains [21,22] to form M -leg ladders. While there are many other possible choices, the SSH chain is one of the simplest systems with non-local, symmetry-protected topological edge states. Furthermore, its properties can be studied experimentally using cold atoms in optical superlattices [23][24][25]. Open boundary conditions can be implemented using an optical box potential [26,27]. As we will show below, the operational entanglement that results from the topological edge states can be observed in this system using modern experimental techniques.
Let J, δ s and z be constants that indicate hopping between sites, and a s j † (b s j † ) creation operators of an 'a' ('b') spinless fermion on chain s, in unit cell j. The unit cell consists of 2M elements, two elements on each of the M chains. The Hamiltonian of the model in second quantization is then given by with J = 1 for the remainder of the paper. A visualization of the SSH ladder for M = 3 with open boundary conditions is shown in Fig. 1.
The non-interacting SSH ladder, a member of the BDI symmetry class, has three non-spatial symmetries. These symmetries are time reversalT , charge-conjugationĈ, and chiralŜ, see Ref. [13]. The relevant symmetry here is the chiral symmetry.Ŝ is defined in terms of its action on annihilation operators asŜ There are also additional non-spatial symmetries present when certain restrictions are placed on the parameters. In particular, additional chiral symmetries enable additional topological invariants [28] which are discussed further in App. C. Topology.-In order to find the number of edge states of the M -leg SSH ladder for a given set of parameters, we have to analyze its topological phase diagram. To define the winding number I for our system with chiral symmetry, we follow Ref. [29].
First, we define a unitary matrix U S for the chiral sym-metryŜ from the conditionŜψ nŜ −1 = m (U S ) * n,m ψ † m , where ψ n are the a s j and b s j operators in an arbitrary basis. U S has the property tr U S = 0 and we have defined the phase such that U 2 S = I. U S can be put into a momentum space, block diagonal form represented by U S (k). Let g(k) = H −1 (k) be the matrix representation of the single particle Green's function corresponding to the Hamiltonian H in Eq. (3). Then the topological invariant is given by [29][30][31] which, for non-interacting systems, is equivalent to the winding number as defined, for example, in Ref. [13].
We prove the equivalence of the invariants in App. A. We note, furthermore, that the related Zak phase for a single SSH chain has been measured experimentally in cold atomic gases [23]. Before analyzing the full topological phase diagram for M = 2 and M = 3 leg ladders numerically, we first note the following important analytical results for the general M -leg case: Suppose that δ 1 = δ 2 = · · · ≡ δ and 0 ≤ |z| cos π M +1 < |δ|. (i) If δ < 0, then I = 0. (ii) If δ > 0, then I = 1 for M odd, and I = 0 for M even.
The proof of this result can be found in App. B. We note that this even-odd effect resembles the celebrated result for spin-1/2 ladders [32].
The topological phase diagram based on the invariant defined in Eq. (5) is presented in Fig. 2 for ladders with M = 2 and M = 3. A similar phase diagram for the M = 2 case using a different parameterization can be found in Ref. [31]. When I = 0 we have no edge states and when I = ±1 we have a single filled edge state. The case which we want to concentrate on in the following is I = 2 where we have two filled edge states.
Entanglement.-The ground state of the SSH ladder is described by a pure state wave function ρ = |Ψ Ψ|. We imagine however, that Alice and Bob have access to only a small number of sites at the edges of the system. Note that in contrast to the often studied case of a bipartition, tracing out the rest of the system leads to a mixed state ρ AB which is our starting point. Now that we have a density matrix ρ AB that only describes the two subsystems we are interested in, we can apply the operational entanglement measure defined in Eq. (1) with E[ρ] being a bipartite measure of entanglement for a mixed state. Regardless of the chosen mixed state measure of entanglement, E op is not easy to compute in general. However, for small dimensions of ρ AB a calculation of its matrix elements using correlation functions is feasible [2]. For the case of two edge states (I = 2) considered here the only projected density matrix which will contribute to E op is ρ 1,1 AB . We will call the two modes on one side of the ladder A 1 and A 2 , and the modes on the other side B 1 and B 2 . Next, we define the projection and analoguously P B which project the ground state onto a (non-normalized) state with a single particle in each subsystem, |Ψ 1,1 = P A P B |Ψ . The matrix elements of the 4 × 4 matrix ρ 1,1 AB can now be computed from correlation functions in the projected state, In systems-such as the SSH ladder considered herewhich are Gaussian, we can use Wick's Theorem to turn the multi-point correlation functions into products of two-point correlation functions. Since we already know from the topological phase diagram, Fig. 2, that we need at least a 3-leg ladder to have two filled edge states, we concentrate on this case in the following. The choice of the few sites controlled by Alice and Bob needs to be based on prior knowledge of where the topological edge states are primarily located. Based on numerical results for strong dimerizations, we choose There are many different entanglement measures one can use to quantify the entanglement of a density matrix [33]. Here we use an additive measure of entanglement, the logarithmic negativity [34,35], where ρ T A AB is the partial transpose with respect to Alice and ||ρ T A AB || is the trace norm of the normalized matrix.
AB ] for the 3-leg ladderusing the same dimerizations δ s as in the topological phase diagrams-are shown in Fig. 3. For δ 2 = −0.75 and δ 2 = −0.25, the regions where E neg [ρ 1,1 AB ] ≈ ln 2 coincide with the regions in Fig. 2 with winding number I = 2. The topology of the system is thus directly tied to the operational entanglement which can be extracted from the system by Alice and Bob. However for δ 2 = 0.25, the region with operational entanglement is much smaller than the region with I = 2. In this case, the edge states are not sufficiently localized. Alice and Bob would need to control more sites to extract all of the operational entanglement present in the two edge states. Virtually identical results are obtained if we use the entanglement of formation instead of the logarithmic negativity, see App. E. From an experimental perspective, the most important result however is that the regions with operational entanglement and winding number I = 2 can be identified by simply measuring the generalized number entropy E n , Eq. (2), on a small number of sites only, see bottom panels in Fig. 3. When the edge states are forming, p(n A = 1, n B = 1) increases leading to a decrease of E n . The number entropy can be measured straightfor- wardly in cold atomic gases by single-site atomic imaging as has recently been demonstrated in Ref. [9].
Bell's theorem.-Moving beyond the indirect observation of the entanglement between the edges by monitoring the number entropy, one of the most fundamental ways to prove that two qubits are entangled is to show that a Bell inequality is broken. Here we will choose the Clauser, Horne, Shimony and Holt (CHSH) [18,20] version of Bell's inequality. Let a, a , b, and b be threedimensional vectors. Let σ A/B be the vector of A/B Pauli matrices. Defining . . . 1,1 = tr(ρ 1,1 AB . . . )/p(1, 1), the CHSH inequality reads To show the breaking of the inequality, we choose the vectors a, a , b and b to be in the x-z plane with a·σ A = cos θ a σ A z + sin θ a σ A x and b · σ B = cos θ b σ B z + sin θ b σ B x . We use the representation of the Pauli operators σ A 2 A 2 and similarly for σ B x,y,z . Results for the 3-leg ladder are shown in Fig. 4 with θ ≡ θ a = θ b /2 = θ b /3 and θ b = 0. For δ 2 = −0.75, corresponding to a winding number I = 2 (see top right panel in Fig. 2), the CHSH inequality in the projected state ρ 1,1 AB is broken. Note that while δ 2 = −0.25 and δ 2 = 0.25 also correspond to I = 2, the edge states are not sufficiently localized in these cases to break the CHSH inequality. Using the Fermi-Dirac distribution, we can also evaluate the correlators in the system at finite temperatures [36]. For the example shown in the right panel of Fig. 4, we see that the CHSH inequality remains broken up to temperatures β = 1/T ≈ 200.
Experimental protocol.-Next, we discuss a possible experimental protocol for showing that Bell's inequality is broken. We can relate elements σ x,z A σ x,z B 1,1 to twoparticle correlators of the full many-body state |ψ . Calculating σ z A σ z B 1,1 then amounts to calculating densitydensity correlations which are experimentally accessible [37]. A function such as σ z A σ x B 1,1 , on the other hand, is more difficult to obtain experimentally because it involves measuring correlators such as While for non-interacting systems measuring the single-particle correlation functions might be possible and is sufficient, we can more generally make use of the matrix operation σ z ⊗ σ where R is a π/2 rotation matrix about the y-axis and I the identity matrix. In the following, we show that by time evolving the full many-body state |ψ , we can implement a rotation operator R on the two-site density matrix ρ 1,1 AB . To do so, we use the time evolution operator exp(−iH t) with H defined in Eq. (3), and λ a constant. We now compare the two-site rotated density matrix (I ⊗ R)ρ 1,1 AB (I ⊗ R † ) with the two-site density matrix σ 1,1 AB (t) obtained from the full time-evolved state |Ψ(t) using the fidelity function for density matrices [38][39][40][41] Since we want to rotate around the y-axis, we set λ = −iκ with κ being a real number. We define F 1 (t) = F (ρ 1,1 AB , σ 1,1 AB ) and F 2 (t) = F ((I ⊗ R)ρ 1,1 AB (I ⊗ R † ), σ 1,1 AB ) and show results for these quantities in Fig. 5. For large κ ∼ 10, F 1 and F 2 oscillate out of phase with a maximum fidelity close to 1. This shows that implementing the coupling (10) allows for an effective rotation of the density matrix.
Conclusions.-We have shown that symmetryprotected topological edge states in a system of itinerant particles can be a resource of spatially separated, nonbipartite, operational entanglement which can be transfered to a quantum register of distinguishable qubits. Two filled edge states are then needed to obtain two entangled qubits. While we have used an explicit construction of such edge states based on coupled SSH chains, the connection established here between the topology of the system and the amount of entanglement which can be extracted from its edge states is general. We have shown, furthermore, that the number entropy measured on a few sites only is an indirect probe of the topological and entanglement properties of a system, which is easily accessible in cold atomic gas experiments. Going one step further, we have also demonstrated that the non-bipartite operational entanglement obtained from the projected ground state of the many-body system is sufficient to break Bell's inequality and presented a protocol to measure these strong quantum correlations experimentally.
The authors acknowledge support by the Natural Sciences and Engineering Research Council (NSERC, Canada) and by the Deutsche Forschungsgemeinschaft (DFG) via Research Unit FOR 2316. K.M. acknowledges support by the Vanier Canada Graduate Scholarships Program. We are grateful for the computing resources and support provided by Compute Canada and Westgrid.

Appendix A: Equivalence of invariants
We will manipulate with g(k) = H −1 (k) into an equivalent form. We first note that U S can be put into block diagonal form in momentum space with the blocks represented by U S (k). The property U 2 S = I for the full matrix implies a similar property of the blocks U 2 S (k) = I 2M where M is the number of legs of the SSH ladder. Also sinceŜ is a nonspatial symmetry, tr U s = 0 implies that tr U s (k) = 0 for the individual blocks as well. Then U 2 S (k) = I 2M and tr U S (k) = 0 imply that we can pick a basis such that The operator conditionŜHŜ −1 = H implies for the momentum blocks that U † S (k)H(k)U S (k) = −H(k). This condition implies that in the same basis as (A.2), Then plugging (A.2) and (A.3) into (A.1), we get Now, we use the polar decomposition D(k) = |D(k)|q(k), where |D(k)| is positive definite and q(k) is unitary. We obtain Now we demonstrate that this method of finding the topological invariant is equivalent to the method based on projection operators [12,13]. The formalism starts by writing the Hamiltonian H(k) in the off block diagonal basis Eq. (A.3). Next, we find the column eigenvectors v n (k) of Eq. (A.3). We define where the sum is over all eigenvectors v n (k) with negative eigenvalues. We also define Q(k) turns out to always be of the form Then the invariant is calculated by plugging q (k) into Eq. (A.5) as q(k). Now, we only need to prove that q (k) = q(k). To do so, first define u n (k) as the normalized eigenvectors of |D(k)|. Then the eigenvectors v n (k) with negative eigenvalues are .
Next, we use Eq. (A.9) as the vectors in Eq. (A.6) to obtain the projection operators P (k) and Q(k). We will also use the fact that the u n (k) vectors form a complete basis. We obtain Hence q (k) = q(k).
where the D(k) blocks are M × M matrices. If we define x(k) = |x(k)|e iθ(k) = (1 − δ) + (1 + δ)e −ik , D(k) can be written as Since the symmetryŜ is a non-spatial chiral symmetry, the symmetry operator in Fourier space can be found by replacing j → k n in Eq. (4). Then one applies the Then v ± s are the eigenvectors of H(k) when z = 0. The eigenvalues are E ± s (k) = ±|x(k)|. Now we transform H(k) into the v +/− s basis. Let V be defined by the unit column vectors v where H V 1,1 (k) and H1, 1 V (k) are M × M tridiagonal matrices defined by Without the phase terms, the above matrix has a well known solution for a free particle with open boundary conditions. By modifying the well known open boundary solution with phase factors, we can find the eigenpairs of H V 2,2 (k) = −H V 1,1 (k) and therefore obtain the eigenvec- 14) The We know that the integral in the M odd case is equal to −2π if and only if θ winds around the origin. That happens when δ > 0. This completes the proof.

Appendix C: Additional Chiral Symmetries and Topological Invariants
Here we give examples of additional chiral symmetries that can be defined under certain conditions. Take, for example, the case of the M = 2 SSH ladder. Under the constraint δ 1 = δ 2 , we can define another chiral symmetry operatorŜ 2 . The transformation is defined aŝ One can compute the topological invariant I 2 corresponding to the chiral symmetryŜ 2 . One can switch j → k n to get the corresponding symmetry in momentum space. Then we can useŜψ nŜ −1 = m (U S ) * n,m ψ † m to get the matrix U S2 (k). For this symmetry, we can work in the same basis as Eq. (A.2) for M = 2. In this basis, Then we can plug this directly into Eq. (5) and solve numerically. We find that the invariant I 2 is zero for −1 < δ 1 = δ 2 < 1 and −1 < z < 1.
As another example, we consider the case of the M = 3 SSH ladder. Suppose we have the constraint δ 1 = δ 3 . Then another chiral transformation isŜ Similarly, in the same basis as Eq. (A.2) for M = 3, we obtain (C.4) Fig. 6 shows a comparison of the invariants I from the main text and I 3 defined here under the condition δ 1 = δ 3 ≡ δ. An 'x' is placed in the regions where the analytic result for the invariant I applies. The magnitude of I 3 indicates the number of topological edge states protected by symmetryŜ 3 . Interestingly, for z = 0.9 shown in the right panels of Fig. 6, there is a region where the invariant I 3 is equal to −2 while the invariant I is equal to zero. In this region, there are two topological edge states protected by the symmetryŜ 3 .

Appendix D: Density Matrix Elements
It is helpful to note and similarly for the B operators. Now we simply list a few important equations. The matrix elements not listed are similar.

Appendix E: Entanglement of Formation
Here we consider the entanglement of formation as an alternative measure of bipartite entanglement of a mixed state. Given a density matrix ρ, we define the entanglement of formation as [42] E F (ρ) = min where E is the pure state von-Neumann entanglement and the minimization is taken over all possible ensembles of the form ρ = i p i |ψ i ψ i |.
One of the problems with the entanglement of formation is the numerical difficulty of minimizing over all possible ensembles. However, there is a method of obtaining the entanglement of formation for any two qubit system [43]. First, define a matrixρ = (σ y ⊗ σ y )ρ * (σ y ⊗ σ y ) where ρ * is the complex conjugate of the matrix ρ. Then define R = √ ρρ √ ρ. and let λ i represent the eigenvalues of R in decreasing order. The eigenvalues λ i are also the square roots of the eigenvalues of ρρ. The concurrence is then C(ρ) = max{0, λ 1 − λ 2 − λ 3 − λ 4 }. Define a variable x = 1 2 + 1 2 1 − C 2 (ρ). The entanglement of formation is then E F = −x ln(x) − (1 − x) ln(1 − x). In Fig. 7, the entanglement of formation is compared to the logarithmic negativity for the M = 3 SSH ladder. We can see that the results are virtually the same.