Effects of Boundary on Orbital Magnetization for a Bilayer System with Different Chern Numbers

The real space formalism of orbital magnetization (OM) is an average of the local OM over some appropriate region of the system. Previous studies prefer a bulk average (i.e., without including boundaries). Based on a bilayer model with an adjustable Chern number at half filling, we numerically investigate the effects from boundaries on the real space expressions of OM. The size convergence processes of its three constituent terms $M_{\mathrm{LC}}$, $M_{\mathrm{IC}}$, $M_{\mathrm{BC}}$ are analysed. The topological term $M_{\mathrm{BC}}$ makes a nonnegligible contribution from boundaries as a manifestation of edge states, especially in the case of nonzero Chern numbers. However, we show that the influence of the boundary on $M_{\mathrm{LC}}$ and $M_{\mathrm{IC}}$ exactly compensates that on $M_{\mathrm{BC}}$. This compensation effect leads to the conclusion that the whole sample average is also a correct algorithm in the thermodynamic limit, which gives the same value as those from the bulk average and the $k$ space formula. This clarification will be beneficial to further studies on orbitronics, as well as the orbital magnetoelectric effects in higher dimensions.