Symmetry protected topological (SPT) phases in free fermion and interacting bosonic systems have been classified, but the physical phenomena of interacting fermionic SPT phases have not been fully explored. Here, employing large-scale quantum Monte Carlo simulation, we investigate the edge physics of a bilayer Kane-Mele-Hubbard model with zigzag ribbon geometry. Our unbiased numerical results show that the fermion edge modes are gapped out by interaction, while the bosonic edge modes remain gapless at the (1 + 1)d boundary, before the bulk quantum phase transition to a topologically trivial phase. Therefore, finite fermion gaps both in the bulk and on the edge, together with the robust gapless bosonic edge modes, prove that our system becomes an emergent bosonic SPT phase at low energy, which is directly observed in an interacting fermion lattice model. PACS numbers: 71.10.Fd, 71.27.+a, 73.43.-f B 2 B , z J J (a) (b) BSPT DMI 0 / z J t = 1 / z J t = 2 / z J t = FIG. 1. (Color online) (a) Illustration of AA-stacked honeycomb ribbon (La 1 = 3, La 2 = 3) with periodic (open) boundary condition along the a1 (a2) direction. a1 = (1, 0) and a2 = (1/2, √ 3/2) are the primitive translation vectors. A1, B1, A2 and B2 are the four sublattices within one unit cell. (b) J-Jz phase diagram of bilayer Kane-Mele-Hubbard model. The bosonic SPT (BSPT, red) and dimer Mott insulator (DMI, blue) phases are separated by a bulk transition. The dashed lines inside BSPT denote the J values, above which one can clearly see the exponential decay of the singleparticle Green's function at the boundary from our finite-size calculations. The relative range of such region becomes wider as Jz increases.term [23, 24, 26]. However, as for the physics on the edge, although the field theory and renormalization group anal-