Crystal magnetic response only skin deep? 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.
 A: This is not true--- you made a mistake with the substitution and variation. As I am sure you know, experimentally, it is false, you can get strong bulk magnetic field penetration into any non-superconducting material.
The substitution you made, replacing $B \times r$ with $A$, is not valid for doing variations with respect to A. The electron operators "r" in $B \times r$ are not the same type of thing as the classical function $A(r)$--- this is a formal error. You can't replace $B \times r$ (the classical function $\times$ an electron position operator) with A (the classical function), because then you are ignoring the fact that r is an operator that acts on electron wavefunctions. This is what comes back to bite you when you take the variation.
You should sort it out physically--- in a superconductor, making an A field produces a coherent bosonic current which cancels the A. In a normal metal, you rearrange the Fermi sea into Landau levels, which are then filled up in opposite senses of current, up to the Fermi energy, so that the net current gives a diamagnetic response. It is easiest to say the error in field theory--- the Fermi field operator which has an expected value is a neutral bilinear $\bar{\psi}\psi$ which counts the number of electrons, not the charged bilinear $\psi\psi$ as in a superconductor, which gives a classical charged field response.
The correct variation to take is using the A in the kinetic term for the many-electron Hamiltonian, replacing $p$ with $p-eA$, and then varying with respect to A. In this case, you get the current from the electron motion. The amount of this current is not proportional to A, unless the electrons are coherently superposed bosons. But the reason you got the boson answer in the Fermi case is that you substituted for the psi field as if they were a classical field, and then you do get the same answer as for a classical charged field, and you get the exponential decay.
