Large-scale simulations of vortex Majorana zero modes in topological crystalline insulators

Topological crystalline insulators are known to support multiple Majorana zero modes (MZMs) at a single vortex, their hybridization is forbidden by a magnetic mirror symmetry $M_T$. Due to the limited energy resolution of scanning tunneling microscopes and the very small energy spacing of trivial bound states, it remains challenging to directly probe and demonstrate the existence of multiple MZMs. In this work, we propose to demonstrate the existence of MZMs by studying the hybridization of multiple MZMs in a symmetry breaking field. The different responses of trivial bound states and MZMs can be inferred from their spatial distribution in the vortex. However, the theoretical simulations are very demanding since it requires an extremely large system in real space. By utilizing the kernel polynomial method, we can efficiently simulate large lattices with over $10^8$ orbitals to compute the local density of states which bridges the gap between theoretical studies based on minimal models and experimental measurements. We show that the spatial distribution of MZMs and trivial vortex bound states indeed differs drastically in tilted magnetic fields. The zero-bias peak elongates when the magnetic field preserves $M_T$, while it splits when $M_T$ is broken, giving rise to an anisotropic magnetic response. Since the bulk of SnTe are metallic, we also study the robustness of MZMs against the bulk states, and clarify when can the MZMs produce a pronounced anisotropic magnetic response.