numerical calculation of haldane model arm chair edge states

Question

hello I am trying to numerical simulate the band structure of the one-dimensional periodic arm chair edge states, I use the pybinding model to construct and simulate the arm chair edge states, because we need to use a 4-site lattice to construct the arm chair edges, I believe this causer some problem with the later calculation of edge bands, according to this, the arm chair edge state should be topologically protected, however my simulation suggests otherwise, when the next nearest hopping is tuned to be below 0.5, the edge state is present, but not only in one direction, . is there anything wrong with my code? Thanks below is my code

# this is pybinding implementation for haldane model with zero mass
# showing edge states of zigzag edge and armchair edge.
import pybinding as pb
import numpy as np
import matplotlib.pyplot as plt
from math import sqrt, pi

pb.pltutils.use_style()
a = 0.5  # [nm] unit cell length
t = 1  # [eV] nearest neighbour hopping
t2 = 0.5 * 1.0j  # [eV] next nearest neighbour hopping

mass_eneygy = [0.2, -0.2]

a_cc = a / sqrt(3)

Gamma = [0, 0]
K1 = [-4 * pi / (3 * sqrt(3) * a_cc), 0]
M = [pi / (3 * a_cc) * sqrt(3), pi / (3 * a_cc)]
K2 = [4 * pi / (3 * sqrt(3) * a_cc), 0]


def monolayer_graphene_4atom(on_site_energy=[-0, 0]):
    t2c = t2.conjugate()
    lat = pb.Lattice(a1=[a, 0], a2=[0, a * sqrt(3)])
    lat.add_sublattices(("A", [0, -a_cc / 2]), ("B", [0, a_cc / 2]))
    lat.add_aliases(("A2", "A", [a / 2, a_cc]), ("B2", "B", [a / 2, 2 * a_cc]))
    lat.register_hopping_energies({"n": t})
    lat.register_hopping_energies({"nnA": t2})
    lat.register_hopping_energies({"nnB": t2c})
    lat.add_hoppings(
        # inside the main cell
        ([0, 0], "A", "B", "n"),
        ([0, 0], "A2", "B2", "n"),
        ([0, 0], "B", "A2", "n"),
        ([0, 0], "A2", "A", "nnA"),
        ([0, 0], "B2", "B", "nnB"),
        #     # between neighboring cells
        ([1, 0], "A2", "B", "n"),
        ([0, 1], "B2", "A", "n"),
        ([-1, -1], "A", "B2", "n"),
        #     # next nearest coupling
        ([1, 0], "A", "A", "nnA"),
        ([0, 1], "A2", "A", "nnA"),
        ([1, 0], "A2", "A2", "nnA"),
        ([1, 1], "A2", "A", "nnA"),
        ([1, 0], "A2", "A", "nnA"),
        ([1, 0], "B2", "B2", "nnB"),
        ([0, 1], "B2", "B", "nnB"),
        ([1, 1], "B2", "B", "nnB"),
        ([1, 0], "B", "B", "nnB"),
        ([1, 0], "B2", "B", "nnB"),
    )
    return lat


graphene_4atom = monolayer_graphene_4atom(mass_eneygy)
arm_chair = pb.Model(
    graphene_4atom, pb.rectangle(7), pb.translational_symmetry(False, True)
)

print(arm_chair.system.num_sites)

plt.subplot(131, title="model")
graphene_4atom.plot(hopping={"draw_only": ["n", "nnA", "nnB"]})


plt.subplot(132, title="edge states")
arm_chair_solver = pb.solver.lapack(arm_chair)
arm_chair_solver.set_wave_vector(pi / a_cc / sqrt(3))
arm_prob = arm_chair_solver.calc_probability(28)
arm_prob.plot()


arm_chair_bands = arm_chair_solver.calc_bands([0, -pi / a/sqrt(3)], [0, pi / a/sqrt(3)])
plt.subplot(133, title="edge states")
arm_chair_bands.plot(point_labels=["$-\pi/a$", "$\pi/a$"])

plt.show()