The Hamiltonian for a single electron in a magnetic field reads $$H=\left(\frac{{\bf p}^{2}}{2m_{e}}+q_{e}\phi\right)+\mu_{B}\left({\bf \hat{L}}+g{\bf \hat{S}}\right)\cdot{\bf B}+\frac{e^{2}}{8m_{e}}\left({\bf B}\times{\bf r}\right)^{2}$$ where symmetric gauge ${\bf A}=\frac{1}{2}{\bf B \times r}$ is used and where the last two term are responsible for paramagnetism and diamagnetism respectively.
The many-body Hamiltonian, weak electron-electron interaction considered, reads $${\cal H} = \sum_{\sigma}\int d{\bf r}\,\left[-\left(\psi_{\sigma}^{\dagger}({\bf r})\frac{\hbar^{2}}{2m^{*}}{\nabla}^{2}\psi_{\sigma}({\bf r})\right)+\boldsymbol{\mu}({\bf r})\cdot({\nabla}\times{\bf A})+\hat{n}_{e}({\bf r})\frac{e^{2}}{2m^*}{\bf A}({\bf r})^{2}\right]$$ where $\boldsymbol\mu({\bf r})=\sum_\sigma \mu_B \psi^\dagger({\bf L}+g{\bf S})\psi $ is local magnetic moment density of electron of total electron angular momentum and define $\mathbf M_\text{para} \equiv \langle \boldsymbol \mu\rangle$
Now that we vary $\bf A$, we should obtain the response current $${\bf j}={\nabla}\times{\mathbf M_\text{para}}+\frac{n_{e}e^{2}}{m^*}{\bf A}$$
It is rather strange to me that (i) the last term is essentially the same as the one in London's equation, which gives a penetration depth $\lambda$, and that (ii) the first term is non-zero only on the crystal surface if $M$ is uniform deep inside, which I think is just the magnetization current we saw in textbooks.
If (i) is true, one may conclude that magnetic response for all crystals (neglect interactions) happens in a thin layer near to the surface (a few nanometer) since magnetic field can only penetrate that far. The length scale that electron density drops to zero is a few angstroms, so we are able to neglect the first term in an appropriate region.