The first term \(H_\upsilon \)  indicates the massless Dirac dispersion, where \(\upsilon \)  denotes the Fermi velocity in monolayer graphene. The second term \(H_\gamma \)  denotes the interlayer nearest-neighbor coupling, so that \(H_\upsilon +H_\gamma \)  describes the Hamiltonian for bare BLG [1]}, [2]}, [3]}. The third term \(H_M\)  represents the proximity-induced magnetism. The fourth term \(H_\Delta \)  represents the staggered sublattice potential, which is zero for bare BLG but can appear when BLG is in proximity with other materials [4]}, [5]}, [6]}, [7]}, [8]}, [9]}, [10]}, [11]}, [12]}, [13]}, [14]}. The fifth term \(H_U\)  indicates the interlayer bias (tunable by gate voltage). Without loss of generality, the last term \({\cal H}^{\prime }={\cal H}_R+{\cal H}_I\)  presents the symmetry-breaking induced perturbations from Rashba-type (\(\lambda _R\) ) and Ising-type (\(\lambda _I\) ) spin-orbit couplings, which are the most common spin-orbit interactions induced in graphene-based interfaces [15]}, [16]}, [17]}, [18]}, [19]}. To describe the dependence of magnetism and staggered sublattice potential on layer and sublattice, we parameterize
\(\Sigma _{\tau ,\sigma }=\left(\begin{array}{cccc}1&0&0&0\\0&\beta &0&0\\0&0&\eta &0\\0&0&0&\zeta \eta \\\end{array}\right),~~\tau _{z,\eta ^{\prime }}=\left(\begin{array}{cc}1&0\\0&\eta ^{\prime }\\\end{array}\right),\) 
