Hofstadter butterfly for square finite lattices $n\times n$ with Zeeman term [1]:

$$ \hat{H} = \sum_{ij, \sigma \sigma^{\prime}} t^{ \sigma \sigma^{\prime}}_{ij} \hat{a}^{\dagger \sigma }_i \hat{a}^{ \sigma^{\prime}}_j - \frac{1}{2} g \mu_B B_z \sum \hat{a}^{\dagger \sigma }_i \sigma^z \hat{a}^{\sigma^{\prime}}_i. $$

Model takes into accoun only nnn hopping $t_{nnn}$ and nn hopping $t_{nn}^{ \sigma \sigma^{\prime}}$ with diagonal $t_{nn}$ and non-diagonal part $t_{nn}^{SOC}$ due to spin-orbit term.

Peierls substitutions is taken into account via phase factor, where $a$ is lattice spacing of square lattice:

$$ t_{ij} \rightarrow t_{ij} e^{\gamma_{ij}}, $$

$$ \gamma_{ij} = -2 \pi i \frac{e}{h} a^2 \frac{1}{2} (x_i + x_j) (y_i - y_j).$$

jupyter-notebook, numpy, matplotlib

2x2 lattice with nearest neighbor hopping term:

2x2 lattice with nearest neighbor hopping term (spin-orbit coupling included):

3x3 lattice with nearest neighbor hopping term:

9x9 lattice with nearest neighbor hopping term: