Units related to chemical potential and orbital magnetization I am studying this paper: Physical Review B 74, 024408 (2006) (arxiv)

Abstract
We derive a multi-band formulation of the orbital magnetization in a normal periodic insulator (i.e., one in which the Chern invariant, or in two dimensions (2D) the Chern number, vanishes). Following the approach used recently to develop the single-band formalism [Thonhauser, Ceresoli, Vanderbilt, and Resta, Phys. Rev. Lett. 95, 137205 (2005)], we work in the Wannier representation and find that the magnetization is comprised of two contributions, an obvious one associated with the internal circulation of bulklike Wannier functions in the interior and an unexpected one arising from net currents carried by Wannier functions near the surface. Unlike the single-band case, where each of these contributions is separately gauge invariant, in the multi-band formulation only the sum of both terms is gauge invariant. Our final expression for the orbital magnetization can be rewritten as a bulk property in terms of Bloch functions, making it simple to implement in modern code packages.

and I'm trying to figure out what units are used in this Figure 8:

throughout the paper they seem to use $\mu_0$, which I believe to be the Bohr Magneton, usually $\mu_B$. I am trying to square their units with equation 49
$$ \frac{dM}{d\mu} = - \frac{C}{2\pi  c}$$
where $C$ is the (integer) chern number, here I $|C| = 1$ and $c$ is the speed of light. This formula should be valid in the gap (white):
$ \mu_0 = 0 \rightarrow \mu_1 = 1$ and $E_0 = 0.35 \rightarrow E_1 = -0.55$
If I assume atomic units I get bohrmagneton and eV I get:
$$\frac{dM}{d\mu} = -0.9 \frac{\mu_B}{\text{eV}} = -0.45 \frac{m_{\text{a.u.}}}{\text{eV}} = -12.24 \frac{m_{\text{a.u.}}}{E_{\text{Hartree}}} $$
compared to
$$-\frac{C}{2\pi c} = -\frac{1}{2\pi \cdot137} = 0.00116171491$$
I tried several units system and I can not get anywhere near the scale of their graphs. I know they are using a model and it is somewhat unit less, but this effects the factors in the equations so there should be some underlying unitsystem involve. What units are they using?
 A: Let's take a closer look! Clearly something seems to be off: recall that in classical electrodynamics the orbital magnetization in the unit cell would be defined as
$$
\mathbf{M} = \frac{1}{2 \Omega} \int_\Omega \text{d}V~\mathbf{r} \times \mathbf{j},
$$
where $\mathbf{j}$ represents the current density (and where $\Omega$ is the unit cell volume). This implies that magnetization has the units of $[\mathbf{M}] \overset{!}{=} \text{Electrical Current}$ in $d=2$ dimensions.
Please note also that the Bohr Magneton $\mu_\text{B}$ is the unit of magnetic moment while $\mathbf{M}$ represents the magnetic moment per volume (surface area). From a direct comparison of the two expressions it is clear that equation (4) from your paper misses the charge $e$ and has an additional factor $c$ (pointing at a different unit system, most-likely the Gaussian system). 
In fact, if you take a look at this paper on the Berry Phase Effects on Electronic Properties, you will get from expression (5.23) the S.I. formula:
$$
\left( 
\frac{\partial M}{\partial \mu}
\right)_{B,T}
= - \frac{e}{\hbar} \frac{\mathcal{C}}{2\pi},
$$
where $\mathcal{C}$ is the first Chern number. Let 
$$[\xi]=\text{Energy}$$
represent the energy scale of your system. If you are now working with a natural unit system, i.e., $e=\hbar=c=\xi=1$ you are effectively expressing your magnetization in units of $e \xi / \hbar$.
As you can see, there will be no arbitrary scale left in this equation and this is why I would deduce that
$$
\mu_0 \equiv \frac{1}{2 \pi}.
$$
This choice is convenient, since the slope is then directly provided by the Chern number, i.e., within the dimensionless system you have
$$ \frac{1}{\mu_0}
\left( 
\frac{\partial M}{\partial \mu}
\right)_{B,T}
= - \mathcal{C}.
$$
Addendum
By checking the image you provide you will find that it takes $\eta\approx 6.76$ times more pixels to represent the number 1 along y than along x. Performing a linear regression:

So you extract the Chern number
$$
\mathcal{C} \approx \frac{6.58}{\eta} = 0.97,
$$
which is what you expect.
