As I wrote in the comments, there is not really a good way of arriving at these conclusions for a general aromatic molecule. The energy differences at play in these systems are so small that qualitative arguments will often be insufficient to get the order right. The diamagnetic term, which can be qualitatively rationalized by arguments based on electron density, for non-H atoms, is dominated by the paramagnetic term, which is, in some sense, an orbital response term due to the application of the magnetic field of NMR. Qualitatively reasoning about the paramagnetic term is possible in some special systems like PhLi or pyridine, but the accuracy of these qualitative arguments is insufficient for general substituted benzenes. The paramagnetic term can be calculated by electronic structure packages, but I assume that is not really the question.
I wrote some lecture notes a few years ago for a lecture I gave on the subject of diamagnetic and paramagnetic terms in NMR. I tried to pare it down a bit but they are still rather long. I will start with the qualitative quasiclassical conclusions and then have the derivation.
Qualitative Explanation of Diamagnetic and Paramagnetic Terms
The diamagnetic and paramagnetic terms in NMR shifts arise from considering how the electron motion is affected by an applied field. The terms diamagnetic and paramagnetic only indicate the sign of the magnetic response, opposed to the field and with the field, respectively (it does not have anything to do with unpaired electrons, which is the place most chemists have heard these words before). The diamagnetic term can be thought of being like Lenz's Law in classical electromagnetism; applying a field induces a current which opposes the field. This is sort of (with excessive handwaving) why electron density correlates with how shielded the nuclei are: with less electron density around, the amount of current is less, the induced field opposing the applied field is less, and the effective field is larger.
For the paramagnetic term, think of an electron moving (remember this a quasiclassical picture, but it's not wrong) in a $p$ orbital. It has angular momentum $l=1$, including an $m_l$ value of $+1$ which would generate a huge field and an $m_l$ value of $-1$ which generates the same field and in the opposite direction (also $0$ but that's not important here). Hence, in the absence of a field there is no field due to this orbital angular momentum being ``quenched'', that is $m_l=\pm 1$ are represented equally in the ground state wavefunction and exactly oppose each other. If a field is applied though, the $m_l$ aligned with the field (let's say $+1$) is stabilized while the $m_l=-1$ is destabilized and so there is more of one than the other in the ground state wavefunction at that angular momentum is unquenched. This results in the addition of a field in the same direction as the applied field, making it larger. This is the origin of the much larger range of shifts for $\ce{C}$ as opposed to $\ce{H}$. Calculating the magnitude of the paramagnetic term is not super straightforward and so qualitative arguments to try to rationalize the order of simple aromatic molecules is hard to do (examples like PhLi or pyridine are possible because there is an obvious paramagnetic perturbation, the lone pair). The amount of paramagnetic shifting tells us something about the low-lying vacant orbitals of the system and so can be a useful tool for predicting reactivity in some cases, but predicting those shifts without calculating them explicitly is not readily possible.2
In hydrogen, which is essentially only an $s$ orbital (which has an angular momentum of 0), there is no paramagnetic term (because there is no orbital angular momentum to unquench) and all of the shifts are due to the diamagnetic term. For heavier elements, like C, where $p$ orbitals get involved, there is a paramagnetic term and it is dominant.
Derivation of the Diamagnetic and Paramagentic Terms
Classical electromagnetism tells us that the magnetic field $\textbf{H}$ generated by a charge $q$ moving with velocity $\mathbf{v}$ at a point $\mathbf{r'}$ away is
$$\textbf{H}=\frac{q}{c}\frac{\mathbf{v}\times\mathbf{r'}}{r'^3}$$
Since we care about the field at the nucleus, we set the coordinate system origin at the nucleus and consider the field at the origin of a charge a distance $\mathbf{r}$ away. Then $\mathbf{r'}=-\mathbf{r}$ and we have
$$\textbf{H}=\frac{q}{c}\frac{\mathbf{r}\times\mathbf{v}}{r^3}=\frac{q}{mc}\frac{\mathbf{r}\times\mathbf{mv}}{r^3}=\frac{q}{mc}\frac{\mathbf{L}}{r^3}$$
We'll need this later. Let's consider though: the above implies that $\textbf{H}$ is on the order of $\approx \beta/r^3$, where $\beta$ is the Bohr magneton and $(1/r^3)$ for a 2p orbital is about 0.25 Å and so $H$ is about 600,000 gauss, which is massive! This should bring up the question: if non-zero $\mathbf{L}$ leads to huge magnetic fields why don't these fields completely destroy our NMR signal? Are there magnetic fields present in diamagnetic molecules when there are no magnetic fields present? I've sort of already hinted at the answer. No, there aren't such huge fields. Loosely, this is because the field generated by the electrons circulating in the $m_L = +1$ level exactly cancels the field of the electrons in the $m_L=-1$ because both are equally represented in the ground state (as they must be without a magnetic field present). This quenching of angular momentum occurs in all cases where the orbitals are real valued, such as if they are solutions of a Hamiltonian in the absence of a magnetic field (so nothing with unpaired spins, which carry a magnetic field). If we turn on an $\textbf{H}_0$ field though, one of the directions of circulation will be favored over the other (the $m_L=-1$) and the wavefunction will adjust so that the ground state has a slight net circulation in the favored direction. In this way, an applied magnetic field unquenches a small amount of the orbital angular momentum. This gives rise to the paramagnetic contribution.
Let's develop this mathematically. First, we should show that our semi-classical picture is a valid way to consider this quantum mechanical system. Fundamentally, a magnetic field $\textbf{H}$ is given by the curl of some vector potential $\textbf{A}$. We ultimately have two fields to consider: $\textbf{H}_0$, the applied field, and $\textbf{H}_n$, the field generated by the nucleus of the atom, each with their associated vector potentials.
\begin{align}
\textbf{H}_0&=\nabla\times \textbf{A}_0\\
\textbf{H}_n&=\nabla\times \textbf{A}_n
\end{align}
We're going to neglect the magnetic fields of the nucleus for this discussion. Obviously, they are important, since the perturbation of the nuclear spin energy levels are mediated by the coupling of the nuclear magnetic fields to the electronic magnetic fields, but we first need to determine what those electronic fields are.
As you may recall from classical electromagnetism, there is more than one vector potential that produces a given field. Specifically, if $\textbf{H}=\nabla\times\textbf{A}$, another vector potential $\textbf{A}' = \textbf{A}+\nabla\phi$, where $\phi$ is a scalar function, gives the same field, since the curl of a gradient of a function is always 0. The transfomation $\textbf{A}\rightarrow\textbf{A'}$ is called a gauge transformation. We need all physical results that we calculate to be gauge-invariant.
So now, how does an applied magnetic field show up in our Hamiltonians. It appears by changing $(\hbar/i)\nabla\rightarrow(\hbar/i)\nabla-(q/c)\textbf{A}$. So the Hamiltonian becomes
$$\mathcal{H} = \frac{1}{2m}\left(\mathbf{p}-\frac{q}{c}\textbf{A}\right)^2+V$$
If we were to calculate expectation values of $\hbar/i\nabla$ with different gauges, we would find that they aren't equal and so $m\mathbf{v}$ is no longer $\hbar/i\nabla$ but is $\hbar/i\nabla-\frac{q}{c}\textbf{A}$. Similarly, $\mathbf{L}$ is now
$$r\times\left(\frac{\hbar}{i}\nabla-\frac{q}{c}\textbf{A}\right)$$
We define something that will be useful to us throughout this discussion, the current density $\mathbf{j}(\mathbf{r})$.
\begin{align}\label{current}\tag 1
\mathbf{j}(\mathbf{r})=\frac{q}{2m}\frac{\hbar}{i}[\psi^*\nabla\psi-\psi\nabla\psi^*]-\frac{q^2}{mc}\textbf{A}\psi^*\psi
\end{align}
Note that it is a function of position and real-valued. It is also $q$ times the quantum mechanical probability current and gauge-invariant. If we assume $\psi$ is a solution of the time-dependent Schrodinger equation, then
$$\text{div} \mathbf{j}+\frac{\partial \rho}{\partial t} = 0$$
where $\rho = q\psi^*\psi$ (that is, the charge distribution). For pure states, $\psi^*\psi$ is time-independent and so $\text{div} \mathbf{j}=0$, and hence behaves as a classical current density. One can do some more math and show that it is possible to express the chemical shift, a quantum mechanical phenomenon, entirely classically as the coupling of the the nuclear magnetic moment to the effective magnetic field at the nucleus due to the circulating density $\mathbf{j}_0$.
\begin{align}\label{chemshiftE}\tag 2
E_\text{pert}=-\mathbf{\mu}\cdot\left[\frac{1}{c}\int{\frac{\mathbf{r}\times\mathbf{j}_0(\mathbf{r})}{r^3}d\tau}\right]
\end{align}
So if we can calculate $\mathbf{j}_0$ we have done the hard part of finding the chemical shift. Let's turn to that now.
Now, the Hamiltonian for our electron is
$$\mathcal{H} = \frac{1}{2m}\left(\mathbf{p}-\frac{q}{c}\textbf{A}_0\right)^2+V$$
We multiply out the binomial and rearrange to get
$$\mathcal{H} = \frac{p^2}{2m}+V-\frac{q}{2mc}(\mathbf{p}\cdot\textbf{A}_0+\textbf{A}_0\cdot\mathbf{p})+\frac{q^2}{2mc^2}A_0^2$$
Let's assume we have the eigenfunctions of $\mathcal{H}_0=\frac{p^2}{2m}+V$ in the absence of the field (i.e. $\mathcal{H}_0\psi_n=E_n\psi_n$). Our goal is to calculate $\mathbf{j}_0$, which has a term linear in $\textbf{A}_0$, which looks like
$$\textbf{A}_0=\frac{1}{2}\textbf{H}_0\times\mathbf{r}$$
That's one choice of gauge, but in any gauge $\textbf{A}_0$ is proportional to $\textbf{H}_0$. So if we want $\mathbf{j}_0$ right to first order, we need $\psi$ at least accurate to first order in $\textbf{H}_0$ (the $\psi^*\psi$ in the last term of $\mathbf{j}_0$ can be the unperturbed wavefunction since that term already has a $\textbf{A}_0$). So let's do some first-order perturbation theory! In case you forgot your first-order PT, we first need to define the perturbing Hamiltonian. We'll use the linear part of our Hamiltonian above.
$$\mathcal{H}_1=-\frac{q}{2mc}(\mathbf{p}\cdot\textbf{A}_0+\textbf{A}_0\cdot\mathbf{p})$$
Our first-order corrected ground-state wavefunction (this is just standard first order PT) is
$$\psi'_0=\psi_0+\sum_n{\frac{\langle n|\mathcal{H}_1|0\rangle}{E_0-E_n}\psi_n}$$
where $E_0$ is the ground state energy, $E_n$ is the energy of the $n$th excited state $\psi_n$. Let's define
$$\epsilon_{n}=\frac{\langle n|\mathcal{H}_1|0\rangle}{E_0-E_n}$$
So that
$$\psi'_0=\psi_0+\sum_n{\epsilon_{n}\psi_n}$$
Substituting this $\psi'_0$ into the equation for $\mathbf{j}$ (equation \ref{current} above), we obtain
$$\mathbf{j}_0(\mathbf{r})=\frac{\hbar q}{2mi}(\psi_0^*\nabla\psi_0-\psi_0\nabla\psi_0^*)+\frac{\hbar q}{2mi}\sum_n\left[{(\psi_0^*\nabla\psi_n-\psi_n\nabla\psi_0^*)\epsilon_n+(\psi_n^*\nabla\psi_0-\psi_0\nabla\psi_n^*)\epsilon_n^*}\right]-\frac{q^2}{mc}\textbf{A}_0\psi^*\psi$$
The first term
$$\frac{\hbar q}{2mi}(\psi_0^*\nabla\psi_0-\psi_0\nabla\psi_0^*)$$
is the current flowing in the absence of a magnetic field. If our angular momentum is quenched as we described above (as in closed-shell molecules), then this is 0. If we take the $\psi_n$ as real, then
$$\mathbf{j}_0(\mathbf{r})=\frac{\hbar q}{2mi}\sum_n\left[{(\epsilon_n-\epsilon_n^*)(\psi_0\nabla\psi_n-\psi_n\nabla\psi_0)}\right]-\frac{q^2}{mc}\textbf{A}_0\psi^*\psi$$
We're gonna see that the first term is the paramagnetic current, related to $\epsilon_n$, which recall is the matrix element coupling the ground state to excited states by the angular momentum operator, and the second term is the diamagnetic current, always present and independent of $\textbf{L}$. Let's go back to that choice of the vector potential $\textbf{A}_0$ we made before
$$\textbf{A}_0=\frac{1}{2}\textbf{H}_0\times\mathbf{r}$$
Subbing that into the perturbing Hamiltonian gives
\begin{align}
\mathcal{H}_1&=-\frac{q}{2mc}(\mathbf{p}\cdot\frac{1}{2}\textbf{H}_0\times\mathbf{r}+\frac{1}{2}\textbf{H}_0\times\mathbf{r}\cdot\mathbf{p})\\
&=-\frac{q}{2mc}(\textbf{H}_0\cdot(\mathbf{r}\times\mathbf{p}))\\
&=-\frac{q}{2mc}(\textbf{H}_0\cdot\mathbf{L})
\end{align}
Where that $\mathbf{L}$ is in the absence of a field. Then we have that
$$\epsilon_{n}=\frac{\langle n|\mathcal{H}_1|0\rangle}{E_0-E_n}=\frac{qH_0}{2mc}\langle n|\mathbf{k}\cdot\mathbf{L}|0\rangle$$
where I assumed that $\textbf{H}_0$ is some magnitude $H_0$ along some unit vector direction $\mathbf{k}$. Now, going back to our chemical shift energy (equation \ref{chemshiftE}), we want the field that this current $\mathbf{j}_0$ generates. This is
$$\textbf{H}=\frac{1}{c}\int{\frac{\mathbf{r}\times\mathbf{j}_0(\mathbf{r})}{r^3}d\tau}$$
Substituting $\mathbf{j}_0$ and $\mathcal{H}_1$ in gives
$$\textbf{H}=H_0\frac{q^2\hbar}{2m^2c^2}\sum\frac{\langle 0|\sum_j\frac{\mathbf{L}_j}{r_j^3}|n\rangle\langle n|\mathbf{k}\cdot\mathbf{L}|0\rangle+\langle 0|\mathbf{k}\cdot\mathbf{L}|n\rangle\langle n|\sum_j\frac{\mathbf{L}_j}{r_j^3}|0\rangle}{E_n-E_0}-\frac{q^2}{2mc^2}H_0\langle 0|\sum_j\mathbf{r}_j\times(\mathbf{k}\times\mathbf{r}_j)|0\rangle$$
where $\mathbf{k}$ is the direction of the applied field and the $j$ subscripts denote electron number in the wavefunction.
For an $s$ orbital, all matrix elements with $L$ are 0 and so the first term is 0. This leaves only the second term which is opposed to the applied field. For anything with higher angular momentum than $s$ orbitals, the first term will be nonzero and aligned with the field. This term arises due to unquenching of the orbital angular momentum in the presence of a magnetic field and hence the direction of the paramagnetic field was determined by the projection of angular momentum which aligns with the field, rather than arising out of Lenz's Law. Said yet a different way, the orbital was already a magnet, it's just that we couldn't see it before. The applied field does not induce the quantum mechanical current flow or paramagnetic field (in the sense of classical EM), merely reveals it by mixing more of the aligned magnet into the ground state wavefunction than the anti-aligned one.
Putting some reasonable numbers into these equations gives that the paramagnetic term is larger at distances less than 1 Å from the nucleus, while the diamagnetic term dominates outside this range. As we mentioned above, the distance important for chemical shifts (average $r^3$ for a $2p$ orbital) is about 0.25 Å and so the paramagnetic term will be larger for NMR related problems. On the other hand, these paramagnetic terms play practically no role in macroscopic magnetic moments of molecules.
References:
- Slichter, Charles P. Principles of Magnetic Resonance (Harper and Row Publishers, 1963).
- Gordon, C. P., Yamamoto, K. , Liao, W., Allouche, F., Andersen, R. A., Copéret, C., Raynaud, C., Eisenstein, O. Metathesis Activity Encoded in the Metallacyclobutane Carbon-13 NMR Chemical Shift Tensors. ACS Cent. Sci., 2017, 3 (7), 759–768. DOI: 10.1021/acscentsci.7b00174
