Time-reversal invariant topological skyrmion phases

Topological phases realized in time-reversal invariant (TRI) systems are foundational to experimental study of the broader canon of topological condensed matter as they do not require exotic magnetic orders for realization. We therefore introduce topological skyrmion phases of matter realized in TRI systems as a foundational step towards experimental realization of topological skyrmion phases. A novel bulk-boundary correspondence hidden from the ten-fold way classification scheme is revealed by the presence of a non-trivial value of a $\mathbb{Z}_2$ spin skyrmion invariant. This quantized topological invariant gives a finer description of the topology in 2D TRI systems as it indicates the presence or absence of robust helical edge states for open boundary conditions, in cases where the $\mathbb{Z}_2$ invariant computed with projectors onto occupied states takes a trivial value. Physically, we show this hidden bulk-boundary correspondence derives from additional spin-momentum-locking of the helical edge states associated with the topological skyrmion phase. ARPES techniques and transport measurements can detect these signatures of topological spin-momentum-locking and helical gapless modes. Our work therefore lays the foundation for experimental study of these phases of matter.

The set of mappings from the full Brillouin zone to the space of an observable O generally exhibits topological sectors if the mappings correspond to a non-trivial homotopy group, however [42,43].The expectation value of O corresponding to such a non-trivial homotopy group can wind over the Brillouin zone to yield quantized, non-trivial topological charge for mappings in topological sectors, this topological charge is, in general, distinct from the topological charge due to winding of the projectors onto occupied states.A physically relevant degree of freedom that has been stud-ied in previous works is spin [42][43][44], though only recent work [42,43] considers cases where the spin topological invariant is not locked in value to the projector topological invariant [44].These topological phases associated with the spin degree of freedom specifically are known as topological skyrmion phases of matter [42], which are now understood to be lattice counterparts and first evidence of quantized transport of magnetic skyrmions, a quantum skyrmion Hall effect [45].Notably, the quantization of topological charge computed as the winding of the ground state spin expectation value over the full BZ is guaranteed whenever the minimum magnitude of the ground state spin expectation value is finite, even when spin is not a conserved quantity [42,43].

II Models
In the following, we consider tight-binding Hamiltonians which preserve time-reversal symmetry (TRS) and realize skyrmions in spin textures over the Brillouin zone.To construct such Hamiltonians, we combine a TRS-breaking Hamiltonian, which realizes a topological skyrmion phase, with its time-reversed partner, and couple these two Hamiltonians with additional spin-orbit coupling (SOC) terms.This approach takes inspiration from construction of some Hamiltonians describing quantum spin Hall insulator phases, by pairing a Chern insulator with its time-reversed partner, and then adding non-negligible SOC terms [6,7,9].A similar construction was used to realize first examples of helical topological skyrmion phases [42] characterized by non-trivial skyrmion numbers computed for mirror subsectors, equal in magnitude and opposite in sign.The present work goes beyond this past study in that block-diagonalization is not required, and time-reversal symmetry is present.While topological skyrmion phases of matter are realized by three-band Bloch Hamiltonians [45], the simplest cases of topological skyrmion phases are, arguably, in four band models with generalized particle-hole symmetry [42,43], we consider an eight band model.The additional degree of freedom has to be given by the spin of the electron since the four band models possess no TR symmetry.We therefore reinterpret the skyrmion for each spin sector as forming in the texture of a pseudo-spin degree of freedom over the Brillouin zone.This is consistent with past work on topological skyrmion phases: ultimately, any pseudo-spin with the appropriate set of deformations to yield a non-trivial homotopy group has the potential to realize a topological skyrmion phase [45].
Explicitly, we consider the following Hamiltonian for the TRI topological skyrmion phase: where s µ , τ µ , σ µ , µ = 0, 1, 2, 3 are Pauli matrices acting on the spin, particle-hole and pseudo-spin degrees of freedom respectively.We have defined ϵ(k) = 2+M −t(cos k x +cos k y ) where t, M are real numbers describing hopping and staggered potential amplitudes.We further restrict ourselves to spin-triplet pairing in the pseudo-spin sector ∆(k) = ∆ 0 (sin k x s τ y σ 0 − sin k y s z τ x σ z ), where ∆ 0 is the pairing strength, and spin orbit coupling term V SOC = cs x τ 0 σ y .
Given the high symmetry of the Hamiltonians, there is some possibility of other symmetries protecting other topological phases besides those characterized by the Z 2 projector invariant of the QSHI, or the skyrmion invariant.We therefore also include the term V bulk , which may contain additional perturbations to break all symmetries save for particle-hole symmetry and time-reversal symmetry as discussed below.For negligible SOC, this Hamiltonian consists of two topological skyrmion phase Hamiltonians, which possess total Chern number C tot of +2 or −2 and skyrmion number Q of −1 or +1 for half-filling, respectively, depending on the spin sector.For non-negligible SOC, this corresponds to a trivial value for the Z 2 projector topological invariant of the quantum spin Hall insulator phase, and non-trivial value for the Z 2 skyrmion topological invariant.
To explore the consequences of the non-trivial skyrmion topology, we additionally consider a counterpart TRS Hamiltonian constructed from TRB topological skyrmion phase Hamiltonians with even total Chern number and even skyrmion number: As these two Hamiltonians each possess PHS, they can be written in the Nambu basis, to describe superconductors at mean-field level using Bogoliubov de Gennes formalism [46,47].Eq. 1 and Eq. 2 are invariant under the particle-hole operators C = s z τ x σ z K and C 2 = τ x K, respectively, which each square to one.By construction, each model is also invariant under the spinful time-reversal operation T = −is y K.

III Topological characterization
To characterize the topology of these models, we employ a generalization of the skyrmion invariant to TRI-systems.The use of the total skyrmion number to characterize the topology fails, since, in analogy to the Chern number, it is always zero when spinful TRS is present, as proved in appendix S1.In the case of the quintessential 2D TR-invariant systems classified by the ten-fold way, one can construct a non-zero topological invariant in multiple equivalent formalisms.One such path involves the usage of Kramer's theorem to express the Hamiltonian in terms of the helicity basis or Kramers pair basis, where the Hamiltonian is block diagonal.This can be done in the case of negligible SOC.Then, for each Hamiltonian one can calculate the Chern number and define a parity index ν = (C − C T R )/2 mod(2) = (−1) C .The parity index is equivalent to the celebrated Z 2 invariant of Kane-Mele [6,9].
In analogy to this parity invariant, we may compute a skyrmion number for each spin sector, characterizing a topological pseudo-spin texture over the Brillouin zone.The explicit form of the spin operators is given in appendix S2 and a schematic visualization of this mapping is shown in Fig. 1.Once these operators are defined, the topological skyrmion invariant can be calculated for each spin sector so as to obtain two Q I , Q II ∈ Z integers.As detailed in appendix S2, these quantities are opposite in sign and equal in magnitude so that it is natural to define a TR-skyrmion invariant or skyrmion parity given by ν Q = (Q I − Q II )/2 mod(2).This invariant, which now fully takes into account the TR-invariant nature of the system, will be shown to be linked to a novel bulk-boundary correspondence when open boundary conditions in x or y are applied.It will be further shown numerically that the helical edge states present when the TR-skyrmion invariant is non-trivial are robust to disorder and local perturbations.Furthermore, the remarkable behaviour displayed by this phase doesn't rely on crystalline point group symmetries, as will be shown by explicitly breaking all symmetries except for time reversal and particle-hole.

IV Results
Starting with the simplest case of negligible spin-orbit coupling, we consider two uncoupled Hamiltonians labeled by spin sector, which each possess four bands, with a welldefined Chern number and skyrmion number as described in Liu et al. [43].For H 1 (k) the total Chern number of each sector takes only two non trivial values C I = ±2 and C II = −C I because of TRS.As the topological invariant for a quantum spin Hall insulator with non-negligible SOC according to the ten-fold way classification scheme is ν = (C I − C II )/2 mod(2), the value of ν for our Hamiltonian is then always trivial as ν = 2 mod(2) = 0 mod 2.
In contrast, Q I = ∓1 in the regions of the phase diagram where C I = ±2.Therefore, the skyrmion number for the TRI system, ν Q , evaluates to 2) for these regions of the phase diagram.The system with non-negligible SOC is therefore in a topologically nontrivial phase according to the skyrmion invariant, even when the projector invariant indicates the system is topologically trivial.
To explore the consequences of the non-trivial skyrmion number paired with trivial projector topological invariant, we consider the case of trivial projector invariant and trivial skyrmion number by comparing results for the Hamiltonian given by Eq. 1 with those for the Hamiltonian given by Eq. 2.
In the case of Eq. 2, C I = ±4 in the topologically non-trivial regions of the phase diagram, corresponding to ν = 0 when SOC is non-negligible, and Q I = ∓2, yielding a trivial ν Q = 0 as well.Comparing these two cases, we observe additional structure that is captured by the skyrmion invariant, which is richer than the ten-fold way classification scheme.
The previous analysis is valid when no SOC term is present.For ν even, it is known that as long as V SOC doesn't break TRS and doesn't close the minimum direct bulk energy gap, then the invariant will always be trivial and the system is expected to be topologically trivial.We therefore explore the effects of non-negligible SOC in systems characterized by the ν invariant, but also the previously-unidentified skyrmion invariant ν Q .
Having defined this TR-skyrmion invariant ν Q , we compute phase diagrams characterizing the skyrmion topology for Hamiltonian Eq. 1, which are shown in Fig. 2. ν Q is shown in Fig. 2 a) for ∆ 0 = 0.5, as a function of SOC constant c and mass parameter M .To understand the stability of the TRI topological skyrmion phase and corresponding quantization of ν Q , we also show the minimum direct bulk energy gap ∆E in Fig. 2 b) and the minimum spin magnitude |⟨S I ⟩| in Fig. 2  c), respectively, each as a function of c and M .For this parameter set, two regimes are distinguished in the case of H 1 (k): the trivial regime of zero skyrmion number, and the region of Quantization of ν Q survives for finite c and is non-trivial for the regime −2 < M < 0 when each of ∆E and |⟨S I ⟩| are finite.For sufficiently large c, |⟨S I ⟩| goes to zero, and ν Q is finite but unquantized, smoothly approaching zero with further increase of c.The pseudospin texture corresponding to one of these unquantized values of ν Q is shown in Fig. 2 d).This situation is analogous to the loss of quantization of the Hall conductivity in a Chern insulator when the Fermi level intersects bands, rather than bands being completely filled or empty.
For this model, ν Q only changes when ∆E = 0 corresponding to a type-I topological phase transition.This corresponds to symmetry-protection of two sets of spin operators, each with 4 × 4 matrix representation, as opposed to a lower-symmetry model with a single set of spin operators with 8 × 8 matrix representations: we have effectively two four-band toy models for topological skyrmion phases still despite finite SOC terms, and it is previously known that fourband models for topological skyrmion phases with generalized particle-hole symmetry C ′ lack sufficient degrees of freedom to realize the type-II topological phase transition [43], although three-band models without C ′ symmetry realize type-II topological phase transitions quite generically [45].
An example of the skyrmion texture with Q I = +1, with M = −1.3,∆ 0 = 0.5,c = 0.8 is shown in Fig. 3 a),b).The time reverse partner with the antiskyrmion Q II = −1 is shown in Fig. 3 c), d) for the same parameter values.
We now include the effect of a symmetry-breaking term additionally to the normal SOC given by: This extra term in the Hamiltonian breaks all spurious crystalline symmetries of the original Hamiltonian in equation (1).Included in those symmetries is inversion so that inversion is broken in the bulk as well as C ′ symmetry.The only remaining symmetries in the model are particle-hole C and time-reversal T symmetry as well as the combined chiral CT .
Fig. 3 shows a winding within each of the two TR sectors (I and II) that remains even with non-negligible SOC.
An important region of the phase diagram is the zero energy and zero pseudo-spin magnitude regime above approximately c = 1.2.In this case because of the gapless spectrum, the skyrmion number destabilizes and Q I , Q II are not well defined just as the prototypical Z 2 invariant.Remarkably, as shown in Fig 2 d), the pseudo-spin still winds in such a manner that almost everywhere in the BZ it forms a skyrmion.The points of zero spin magnitude seem to not deform the global structure that give rise to the skyrmion texture, although the TR skyrmion invariant is destabilized in this case.It is worth noting that we can also do the same analysis for the second Hamiltonian H 2 (k), in this case the pseudo-spin textures for each TR subsector result in a quantized skyrmion invariant of ±2 and thus not only a trivial projector invariant but also a trivial skyrmion parity.which is shown as a function of SOC in Fig. 4 along with an example of a skyrmion texture.It is worth mentioning that even if the skyrmion parity is trivial across the phase-diagram we still observe quantization of the skyrmion invariant as well as the non-trivial relation Q I = −Q II holding for all parameter values.

V Bulk-Boundary Correspondence
The two Hamiltonians, Eq. ( 1) and Eq. ( 2), each possess a C ′invariant Chern insulator with even Chern number and its TR partner.For Eq. ( 1), the magnitude of this Chern number is 2, and for Eq. ( 2), it is 4.For non-negligible SOC or other perturbations breaking S z conservation, the Z 2 projector invariant calculated from ν = (C−C T R )/2 mod(2) will always be trivial for these Chern numbers even as long as the bulk gap does not close according to the ten-fold way [32,39].One therefore expects the edge states present for zero SOC associated with the Chern insulators will hybridize and gap out immediately for finite coupling between the sectors related by TRS.
More generally, all even Chern number Hamiltonians coupled to their TR partners are predicted to yield a topologicallytrivial phase due to Z 2 classification in these cases.We observe a counterexample to this statement, that all Hamiltonians with even total Chern number coupled to their TR partners are predicted to yield a topologically-trivial phase, for ν Q odd in value, however.We first consider the Hamiltonians (Eq.( 1) and Eq. ( 2) for λ = 0, c = 0 corresponding to negligible SOC and negligible bulk perturbation, respectively.Each exhibits localized edge states in its slab spectrum, that connect bulk conduction and valence bands.
We observe an even number of Kramers pairs at the boundary in agreement with TR symmetry.For a QSHI, finite c corresponding to non-negligible SOC term is expected to gap out the helical edge modes as the topological classification associated with mappings to projectors onto occupied sites reduces from Z to Z 2 .However, in the model considered here, the boundary states appearing from applying OBC to the Hamiltonian (1) do not gap out for finite c and non-negligible SOC.The helical modes instead remain gapless and localized on the edges.Furthermore, the edge states of Eq.( 2) gap out for finite c corresponding to non-negligible SOC, with the energy gap being proportional to c, as expected for a trivial phase.This shows there is some correspondence between additional robustness of the helical edge modes and ν Q odd in value.
We further investigate the mechanism for this unexpected robustness of the helical edge modes for trivial Z 2 projector invariant and non-trivial Z 2 skyrmion invariant by considering λ finite, where λ is the free parameter of a bulk perturbation that breaks all symmetries of Eq. ( 1), except for TR and particle-hole.The results of this perturbation are summarized in Fig. 5. Fig. 5 a) and b) show that, even with all crystalline symmetries broken, the gapless edge states are present.This shows that the gapless edge states are not protected by crystalline point group symmetries of the bulk.
We investigate the robustness of the gapless modes against edge perturbations.Fig. 5 c) and d) show the slab spectra for an additional edge perturbation in Eq. ( 1), which has previously been shown to gap edge states of the QSHI [48]: This perturbation can also be applied for the case of OBC in the x direction by interchanging the x and y labels in its definition.As is clear from the plots in Fig. 5 a)-d) the gapless points may shift, but they remain present even in the present of spin-orbit coupling.
We also investigate the effects of on-site disorder on the edge dispersion and localization of the edge states protected by the skyrmion topology.We add on-site disorder of the form V j s 0 τ z σ 0 for OBC in y, with j being the layer index.This form of disorder preserves time-reversal and particlehole symmetries.The effects on the energy spectrum are shown in Fig. 5 e) and f), respectively.We find that the gapless edge states persist in the presence of a random disorder realization, as shown in Fig. 5 e).To compute corresponding disorder-averaged results, we compute the spectrum of Eq. ( 1) in this slab geometry for each k x , and then combine these en-ergy eigenvalues for each k x into a single larger array, which we then sort in energy.This is computed for each of 100 disorder realizations, and the average of these sorted spectra is shown in Fig. 5 f).We find the edge states of this disorderaveraged spectrum remain gapless, corresponding to each disorder realization simply shifting the crossing point of the helical modes in k x .We now also explore how the non-trivial skyrmion invariant ν Q alters the character of the edge states to prevent hybridisation.Given related work introducing the observableenriched partial trace and characterizing the additional spin-momentum-locking of edge states due to non-trivial skyrmion number [? ], we compute the pseudo-spin expectation values for the edge states for a given TR sector in a slab geometry with either x or y directions open so as to have k y or k x as a good quantum number.We find that, even though non-neglible SOC is present, the edge state pseudo-spin texture has a very distinct pattern depending on the orientation of the edge.For open boundary conditions in the y direction, near the crossing of the edge modes the pseudo-spin S I y , S I z seems to average to zero for each edge, while S I x does not as seen in Fig. 6.Instead the value of S I x on each edge seems to be minus the one of the other edge.
We additionally confirm the nontrivial nature of the edge states by computing the σ xx conductivity as a function of energy as shown in Fig. 6 e).Here the conduction band of the system starts around 0.6 in units of the hopping, since this was set to be one by default, and so a trivial zero conductance would be expected below it.Instead a quantized longitudinal conductivity of 4e 2 /h appears, even in the presence of disorder.Such a signature implies the existance of transmision channels present in the system which necessarily come from edge state transmission since the bulk material is an insulator [49][50][51].Specifically for the case of the QSHI it has been shown that a quantized transport robust to disorder, as in Fig. 6 e), is indeed linked to the presence of topologically protected edge states which for the QSHI carry a 2e 2 /h conductance [52].We indeed observe a quantization to 4e 2 /h consistent with the four edge states per edge which appear when opening boundary conditions as discussed in the previous paragraphs.The robustness to disorder for this ν = 0, ν Q = 1 case is consistent with quenching of SOC by the topological spin polarization of the edge states that is part of the bulk-boundary correspondence of topological skyrmion phases of matter.We may also understand it from the perspective of the lower-symmetry realizations of topological skyrmion phases in study of the quantum skyrmion Hall effect [45]: three-band tight-binding models for topological skyrmion phases of matter without C ′ symmetry realize generalized Thouless pumps, in which spin angular momentum is pumped, rather than electric charge.The spin polarization of edge states observed here corresponds to a C ′ -symmetric and TR-symmetric Thouless pump of spin angular momentum, and gaplessness of edge states is required for this pumping in correspondence with ν Q non-trivial.
However, the most robust astonishing feature associated to the skyrmion invariant is seen upon tracing out the particlehole degree of freedom in each spin sector, by performing the observable-enriched partial trace defined in Appendix S3.That is, we characterize the skyrmion topology by tracing out the particle-hole degree of freedom-in addition to a virtual cut in real-space-to compute a reduced two-point function and corresponding entanglement spectrum computed from the ground state as seen in Fig. 6 f) .In correspondence with the non-trivial skyrmion invariant, this observable-enriched en-tanglement spectrum exhibits helical edge modes.We may therefore understand the spin polarization of edge states in the full system as consequences of the bulk-boundary correspondence of a potentially open subsystem, realized by tracing out the particle-hole degree of freedom.In addition to the previous characterization of the edge states, we finally also explore the remarkable robustness of the gapless helical edge modes for trivial projector invariant and non-trivial skyrmion invariant, by enumerating all possible spin-orbit coupling terms that preserve time-reversal sym-metry and particle-hole symmetry.Out of all possible six momentum-even terms, the gapless helical edge modes are robust against three of these terms, and the other three gap them out.These terms are provided in appendix S4.We distinguish the three SOC terms preserving the gapless helical edge modes from those which gap them out by additional symmetry-protection.
First, we find that certain 4 × 4 submatrices of the Hamiltonian must preserve (CT ) p symmetry (defined in appendix S4), otherwise an effective Rashba SOC term is introduced, permitting finite y component of spin in sector I/II for the edge states in the vicinity of the gapless point(s) in the edge spectrum.That is, the spin-momentum-locking of the edge states is relaxed.As discussed in appendix S4, clockwise spin rotation about the z-axis within sector I/II for one state at the edge is compensated by counter-clockwise rotation of the second state at that edge within the same sector (the C ′ partner of the first state), such that the net spin of these edge states in combination remains polarized in the x-direction.We may compute the observable-enriched entanglement spectrum as defined in Appendix S3 to examine bulk-boundary correspondence in our system upon tracing out the particle-hole degree of freedom of each spin sector to further characterize this edge state spin texture, shown in Appendix S4.We find the spin textures are consistent with gapless states of the spin subsystem propagating along the edge deep within the bulk gap of the entanglement spectrum.Indeed, we examine the edge spectrum of the spin subsystem explicitly to show the nontrivial spin topology is preserved in this lower-symmetry case, but the observable-enriched entanglement spectrum exhibits additional k x -dependent Rashba splitting.This finite y spin component permits hybridisation between edge states previously forced to zero by the stricter spin-momentum-locking for (CT ) p preserved, gapping out the helical edge states protected by the spin topology.
Second, in this system, we have three two-fold degrees of freedom, corresponding to the s, τ , and σ Pauli matrices, respectively.The symmetry operator C ′ , which acts on the τ and σ degrees of freedom only ( acting simply as a constant in the s sector), is required to define two skyrmion numbers, Q I and Q II , computed as winding of spin S I and S II in terms of the matrices discussed in appendix S2, respectively.The precise statements of these symmetry requirements are provided in appendix S4.This additional symmetry protection reflects the fact that the helical edge modes are being protected by topology of subsets of the degrees of freedom, and indicates the concept of symmetry-protection for topological phases must be generalized in light of the topological skyrmion phases of matter.

VI Conclusions
In this paper, we introduce time-reversal-invariant topological skyrmion phases of matter.These are topological phases realized in systems with time-reversal symmetry as required to protect the quantum spin Hall insulator phase, but which are distinct from the QSHI phase.The QSHI is a topological phase of the ten-fold way classification scheme, resulting from mappings from the full Brillouin zone to the space of projectors onto occupied states, while the TRI skyrmion phase arises from topologically non-trivial mappings from the full Brillouin zone to the space of ground state spin expectation values.The TRI skyrmion phases are therefore characterized in the bulk by formation of a skyrmion in the texture of the ground-state spin expectation value over the Brillouin zone in general, rather than in the texture of the projector onto occupied states over the Brillouin zone.
We construct toy models for the TRI skyrmion phase by combining the Bloch Hamiltonian for a C ′ -invariant but time-reversal symmetry-breaking topological skyrmion phase with its time-reversed partner into a Bloch Hamiltonian matrix representation with time-reversal symmetry.We may therefore characterize the TRS skyrmion phase with two skyrmion numbers related to one another by TRS in simpler cases corresponding to a helical topological skyrmion phase.More generally, we may also included terms V SOC coupling these two sectors, in which case the phase is characterized by a single topological invariant ν Q , which is equal to the skyrmion number for one sector, modulo 2. This corresponds to a more general Z 2 classification of ν Q .
Performing an observable-enriched partial trace on the density matrix of the occupied states for a system with open boundary conditions corresponding to a slab geometry to compute the observable-enriched reduced density matrix of the spin subsystem, we find odd ν Q corresponds to topologically-protected gapless modes in the observable-enriched entanglement spectrum.We show this bulk-boundary correspondence of the spin subsystem has consequences for the bulk-boundary correspondence of the full system, demanding robust helical gapless edge states at the boundary of the system for open boundary conditions, even when the system is not in a quantum spin Hall insulator phase according to the Z 2 projector invariant, ν.That is, for ν even and trivial, but ν Q odd and non-trivial, we find topologically-robust, helical gapless edge modes protected by non-trivial ν Q .When each of ν and ν Q is even, topologicallyprotected helical edge modes are absent.
We find the persistence of the helical modes for ν Q non-trivial-when ν is trivial-is related to the strict spinmomentum-locking enforced on edge states in the full system by ν Q non-trivial.This spin-momentum-locking derives from the requirement of gapless helical boundary modes in the observable-enriched slab entanglement spectrum of the spin subsystem due to non-trivial ν Q in the bulk.We add terms to the bulk Hamiltonian such that projectors onto the occupied states possess only time-reversal symmetry T and particle-hole symmetry C. We find that the helical gapless modes observed for ν trivial and ν Q non-trivial are robust for non-negligible V SOC respecting a C ′ symmetry of a subset of the degrees of freedom, as well as lacking a generalized Rashba SOC, which softens spin-momentum-locking at the edge while still yielding the gapless helical edge modes of the spin subsystem after observable-enriched partial trace.Fu-ture work will investigate the potential of TRI six-band models constructed from the lower-symmetry three-band models for topological skyrmion phases [45] to realize topologicallyrobust gapless boundary modes, with the expectation that absence of C ′ symmetry in these models will yield more robust helical edge modes than those observed in the C ′ symmetric cases discussed here, as well as type-II topological phase transitions [42,45].
Notably, these helical gapless modes affect computation of the projector topological invariant: spectral flow winds nontrivially in correspondence with presence of the helical gapless modes, contradicting Z 2 classification of ν.This indicates the projector invariant ν must be generalized to incorporate changes due to ν Q , rather than treated as independent of ν Q .This is also consistent with results of three-band models for topological skyrmion phases, where non-trivial total Chern number as determined by bulk spectral flow does not necessarily yield chiral transport on the edge.These failures relate to reliance of much topological characterization on the flatband limit assumption and topological stability up to closing of a charge gap, which does not hold in general [42,45].
Such topologically-protected gapless helical modes and additional spin-momentum-locking may be observed in transport measurements or spin-ARPES and serve as important signatures and consequences of topological skyrmion phases of matter.The realization of these signatures in time-reversal invariant systems is also highly-desirable for experiment, as the non-trivial topology is realized without the need for particular magnetic orders which generally complicate experimental realization.Our work is therefore significant in bringing study of topological skyrmion phases-including their role in extending topological classification schemes-much closer to experimental realization.

S2 TR skyrmion invariant
In this appendix we focus on defining an invariant that characterizes the topology of the TR skyrmion phase given that the total skyrmion number vanishes.To do that we start with the total spin matrix and consider first the system without TR sector coupling.Because of the particle-hole symmetry this operators take the form: Consider the projector onto the ground state for a given crystal momentum k in the form of: where n = 1, . . .n occ for n o cc the number of occupied bands.Let he time reversal operator be given by T = QK where K is complex conjugation and It is alsop true that for time reversal systems [53]: , where Einstein summation is assumed for α, β which run over all components of the vector |u n (k)⟩ and V nm k = ⟨u n (−k)|T |u m (k)⟩ is a sewing matrix between TR partners.It is then easy to prove how the projector unto the Ground state changes under time-reversal: ,where we used the fact that n V nm k u β n (−k) = Q βα ūα n (k) so that wit some rearranging we obtain: So that we obtain how the projector onto the ground state transforms under time reversal: We are interested not on the projector alone but on the pseudo spin texture which may generically be obtained from the expectation value of some spin-like operators S II µ , S I µ where µ runs over space indices µ = x, y, z, in our case the pseudo-spin operators may be generically given by S I µ = Sµ ⊗ (1 + s z )/2 and S II µ = Sµ ⊗ (1 − s z )/2 .Here we suppose that the basis is such that there is only matrix elements of S I,II µ in one spin sector, where small s µ are Pauli matrices representing the spin degrees of freedom.The matrix Sµ is in principle arbitrary but will turn out to be constrained if TR is satisfied for the textures.We can obtain the pseudo-spin texture by computing the trace: where we assumed the unitary part of the TR operator to act only on the spin degree of freedom via Q = −is y ⊗ I Then considering that the winding number for each pseudo-spin texture can be computed from the effective Hamiltonians: we then have the following relation: This textures are linked then to the pseudo-spin expectation value in a simple way.By using the fact that ρ ↑ (k) = d 0 (k)σ 0 + d ν (k)σ ν (Einstein summation implied) and ρ ↓ (k) are 2 by 2 hermitian matrices generated by Pauli matrices one finds: It is clear then that the pseudo-spin topology is equivalent to the topology of the reduced density matrix for the pseudo-spin subsector defined above, with the final relation: This relationship can be tested numerically and we found an excellent agreement between both textures, an example texture is shown in Fig.
. Additionally the link discussed in the main text between the reduced entanglement spectrum or reduced 2-point correlator is clear from the following derivation.We start with the known non-interacting result for the 2-point correlator which states that the correlation matrix is given by [54]: where ñ∈occ P αβ * ñ (k) gives the matrix elements for the ground state projector P GS (k) of the fully periodic boundary Hamiltonian and we have expressed the real space y correlation matrix as G αβ nm (k x ).Here c † kxnα creates an electron with momentum k x in the n lattice site and in state α.Where α is the local Hilbert space in each site which for us consists of particle-hole,spin and pseudo-spin.Next let us multiply by the unitary considered above and project to the up spin sector so that one obtains the fourier transform of P ↑ U P GS (k)U † .Then proceeding to take the trace over all other degrees of freedom except pseudo-spin σ we have: Furthermore let us trace over N y − L lattice sites, subsystem B, so that we have re reduced correlation matrix which will be a Toepplitz matrix of dimension d loc * L, subsytem A, which we denote GA ↑ (k x ).Now because there is no periodicity in L < N y we have an open boundary condition matrix which we now wish to diagonalize.Since the elements of this matrix are all given by the real space transform of ρ ↑ (k) σσ ′ with no coupling between site L and site 1, for L large enough so that no nearest neighbour hopping exists, we have just the open boundary of the matrix d ↑ 0 (k)σ 0 + 1 2 S I µ k σ µ .Since the Chern number of ρ ↑ (k) is just the Skyrmion invariant and applying the ten-fold way bulk-boundary correspondence for ρ ↑ (k) we conclude that the spectrum of GA ↑ (k x ) will contain Q ↑ chiral gapless states which cannot be gapped out by any perturbation.The same reasoning leads to Q ↓ chiral gapless states for the down projected reduced correlation matrix.
To visualize the Z 2 nature of the invariant we can repeat the same argument but now not taking the projector and tracing out only the particle-hole d.o.f.we obtain the same relation but with Tr τ {U P GS (k)U † } which is now a 4 by 4 matrix in spin and pseudo-spin d.o.f as such the upper 2 by 2 component is just the matrix ⟨↑, σ| Tr τ {U P GS (k)U † } |↑, σ⟩ and the down-spin matrices.These components are related by TR and tied to the Skyrmion invariant since the following relation holds: similar relation holds for the down component and thus both matrices are related by TR symmetry as proven above for the reduced density matrices of the bulk.Moreover, the invariant for each sector Q ↑ = −Q ↓ is then proven to be linked to the topology of Tr τ {U P GS (k)U † } since the Z 2 invariant in this case is just ν The robustness of the reduced entanglement spectrum or equivalently of the reduced correlator spectrum which is shown in the main text against all TR and PH symmetry allowed perturbations indicates that the subsystem only needs the previous symmetries in the full system to maintain TR symmetry.This is clear if the matrix U is the identity since then Tr τ ρ(k) will also have TR symmetry just by using equation (S15) and that the trace does not involve the W matrix so it can get out of the partial trace.The case of an observable which includes a different unitary matrix needs to be treated more delicate since the unitary matrix will mix and entangle degrees of freedom involved in the partial trace.Specializing to our case with U = (1 − σ y )τ z /2 + (1 + σ y )τ 0 /2 we see that the trace now becomes: ,where ρ σ ′ s ′ τ σsτ are the matrix elements of the traced out ground state projector with momentum dependence being implicit, and A τ σ,σ ′ = (τ δ σ ′ ,−1 + δ σ ′ ,+1 )(τ δ σ,−1 + δ σ,+1 ) come from the form of the unitary matrix U considered before in the basis of σ y eigenvectors.Now requiring time reversal symmetry for this reduced ground state projector implies a fix form and relation between the matrix elements for different momentum values.The diagonal terms for example need to fulfill the usual relationship of a TR hamiltonian meaning: this is clearly fullfilled by the TR invariance of the full ground state projector if σ = σ ′ since then the A tensor is the identity.For the antidiagonal case not much can be said unless there is a relationship involving τ ρ(k) which indicates that a particle hole symmetry of the full system is needed since PH always results in an additional negative sign and a reversal of spin.Finally the same reasoning works for the antidiagonal components since TR symmetry implies: ,this can be fulfilled trivially for the σ = σ ′ components by using the TR symmetry of the full system and for the off diagonal in σ the relation will be fulfilled with the additional PH symmetry of the full system which changes the τ sign coming from the A tensor.

S4 Bulk-boundary correspondence and symmetry protection
As discussed in the main text, the edge states encountered for a non-trivial skyrmion parity remain gapless even in the presence of non-negligible spin-orbit coupling and a bulk perturbation breaking crystalline point group symmetries.Here, we investigate in detail the symmetry protection of the edge states and the role of spin-orbit coupling.As a starting point we impose that the system is invariant under T = −is y K and C = s z τ x σ z K.The allowed momentum even spin-orbit terms that couple both up and down spin are then reduced to just six terms: FIG. 1. Pseudo-spin textures of the TR skyrmion phase analyzed in the main text.The arrows represent the pseudo-spin expectation value S I µ k for one TR sector denoted by I, this skyrmion texture has a TR partner (lower skyrmion texture), which is given by the pseudo-spin expectation value S II µ k of the II TR sector.The color is proportional to the polar angle of the texture.The skyrmion numbers correspond to QI = +1, QII = −1 for the model of (1) with M = −1, ∆0 = 0.5,λ = 0.3,c = 0.5.
FIG. 2. a) Phase diagram of the TR-skyrmion invariant νQ for a fix value of ∆0 = 0.5 and bulk perturbation λ = 0.3 of model (1).b)Minimum direct gap as a function of c and M for the same parameter values c) Minimum pseudo-spin magnitude for the first TR or spin up sector.d) skyrmion texture for vanishing magnitude c = 1.5,M = −1,for the first TR sector, the arrows represent the vector (⟨S x I ⟩, ⟨S y I ⟩) while the color represents ⟨S z I ⟩
FIG. 4. a) Skyrmion numbers QI (red) and QII (blue) computed for M = 1, ∆0 = 0.5 as a function of SOC c for the model Hamiltonian of equation (2).b) Minimum direct gap between bands enclosing zero Fermi energy.c), d) Normalized ground state pseudo-spin expectation value plots with t = 1, ∆0 = 0.5,λ = 0.4,c = 0.5.Subfigure c) shows the x − y pseudo-spin texture |⟨(S I x (k), S I y (k))⟩|, while Subfigures d) shows the normalized z component of the pseudo-spin texture ⟨S z I ⟩ FIG. 5. a) Energy spectrum for the model of Eq. (1) with OBC in y and periodic in x with parameters M = −1, ∆0 = 0.5,c = 0.4 length Ny = 200 and λ = 0.0.b) Energy spectrum for the same conditions except λ = 0.3.c) Edge spectrum with OBC in y for M = −1, ∆0 = 0.5,c = 0.5,λ = 0.3 with the edge perturbation discussed in the text V1 = 0.1.d)Same parameters and perturbation but now with OBC in x. e) Edge spectrum with OBC in y for the same parameters as c) but with some random disorder realization respecting particle-hole and TR.d) Same parameters and perturbation but now with OBC in x. f) Energy spectrum as a function of eigenvalue index i for a disorder average over 100 on-site potentials Vjs0τzσ0 distributed with zero mean and standard deviation 0.1.For a system with open boundary conditions in y, periodic in x and parameters M = −1,∆0 = 0.5,c = 0.5, λ = 0.3.
FIG. 6. a) x component of the pseudo-spin expectation value for the edge states of the Hamiltonian in equation 1 with parameters M = −1,∆0 = 0.5, λ = 0.3 , c = 0.5 and open boundary conditions on y.The color is proportional to the localization on either edge computed by taking the probability density at each edge and subtracting them.b) and c) show the y, z components of the pseudospin expectation value for the edge states respectively.d) Probability density for the edge states of the Hamiltonian in equation 1.e) Conductivity along xx as a function of energy for the Hamiltonian in equation 1 with parameters M = −1,∆0 = 0.5, λ = 0.0 , c = 0.5 and full open boundary conditions, random gaussian disorder along the whole sample with standard deviation 0.01 is present.f)2-point reduced correlation spectrum for the Hamiltonian in equation 1 with parameters M = −1,∆0 = 0.5, λ = 0.3 , c = 0.5 and all T , C symmetry allowed perturbations with amplitude 0.5 turned on, for a cut along the y axis and traced over L = 200 lattice sites and particlehole d.o.f.
FIG. S8. a) Energy spectrum for the model of Eq. (1) with OBC in y and periodic in x with parameters M = −1, ∆0 = 0.5,c = 0.5 length Ny = 200 and λ = 0.3 with the edge perturbation discussed in the text V3 of strength 0.2.b) 2-point reduced correlation spectrum for the Hamiltonian in equation 1 with parameters M = −1,∆0 = 0.5, λ = 0.3 , c = 0.5 and perturbation V3 of strength 0.2 turned on, for a cut along the y axis and traced over L = 200 lattice sites and particle-hole d.o.f.c) x component of the pseudo-spin expectation value for the edge states of the Hamiltonian in equation 1 with parameters M = −1,∆0 = 0.5, λ = 0.3 , c = 0.5 and open boundary conditions on y.The perturbation V3 was added with amplitude 0.2.The color is proportional to the localization on either edge computed by taking the probability density at each edge and subtracting them.d) and e) show the y, z components of the pseudo-spin expectation value for the edge states respectively.
Since we still have a non-diagonal SOC matrix element ⟨↑, σ| Tr τ {U P GS (k)U † } |↓, σ ′ ⟩ one would need to calculate the Kane-Mele index for the 4 by 4 matrix Tr τ {U P GS (k)U † } if only time-reversal symmetry is conserved in the subsystem.In this sense applying the reasoning of the previous chapter for the correlation matrix we can see that the helical edge states of opening boundary conditions for Tr τ {U P GS (k)U † } results in a gapless reduced correlation spectrum from Gν (k x ) Tr τ {U G nm (k x ) U † } σs;σ ′ s ′ protected solely by time-reversal symmetry.