I think it helps to analyze this problem in the following two limits: (i) when $E_C \gg E_J$ (the charge qubit limit) and (ii) when $E_J \gg E_C$ (the transmon limit). Since the charge qubit came first, let's begin our analysis there.
If $E_J \rightarrow 0$, the eigenergies are given by $$E_n = 4E_C (n-n_g)^2,$$ where $n$ is the integer number of cooper pairs on the capacitor plate. Note that the energy basis is the (cooper pair) number basis. In other words, the $n^{th}$ excited corresponds to a state with $n$ cooper pairs on the island. By plotting this expression for $E_n$ versus $n_g$ for $n \in {\ldots,~-2,~-1,~0,~+1,~+2, \ldots }$, one obtains Fig. (a) from the <PERSON> paper (apart from the small splittings $\sim E_J$) as depicted below:
But where do the splittings in Fig. (a) come from and why does the magnitude of the splitting appear to decrease as $n$ grows? Let's confine ourselves to the $01$-subspace with $E_J\neq0$ and compute the charge dispersion of the ground state. 976
We find that our Hamiltonian is given by $$ H_{01}(n_g) = \begin{pmatrix} 4E_c(0-n_g)^2 & -E_J/2\ -E_J/2 & 4 E_C(1-n_g)^2 \end{pmatrix}. $$ At $n_g=0.5$, the diagonal entries of this matrix are equal, in other words $$H_{01}(n_g=0.5) = -\frac{E_J}{2}\hat{\sigma}x + E_c \times \hat{1}.$$ The eigenstates and eigenenergies of this matrix are given by: $$|1_q\rangle = (|0\rangle - |1\rangle)/\sqrt{2} \leftrightarrow +E_J/2 + E_C$$ and $$|0_q\rangle = (|0\rangle + |1\rangle)/\sqrt{2} \leftrightarrow -E_J/2 + E_C.$$ At $n_g=0$, we find that $$ H{01}(n_g = 0) = \begin{pmatrix} 0 & -E_J/2\ -E_J/2 & 4 E_C \end{pmatrix}, $$ which has eigenstates and eigenenergies given by: $$|1_q\rangle \simeq |1\rangle \leftrightarrow 4E_C - \frac{E_J^2}{16 E_C}$$ and $$|0_q\rangle \simeq |0\rangle \leftrightarrow \frac{E_J^2}{16 E_C}.$$ Note that we have ignored the effect of the off-diagonal coupling on the eigenstates, but we have included the energy level shifts as they will be useful later. Putting this all together, the ground state charge dispersion $$ \epsilon_0 = |E_0(n_g=0.5)-E_0(n_g=0)| = E_C - E_J/2 \approx E_C. $$ Graphically, $\epsilon_0$ roughly corresponds to the energy at which two adjacent parabolas intersect as measured from 0.
One can apply a similar line of reasoning to compute the charge dispersion of higher excited states. I won't lay out the details, but a key difference is that the diagonal entries for each of the relevant subspaces can no longer be made equal (as was the case for $H_{01}(n_g=0.5)$. Because of this, the energy splittings are reduced to small shifts $\sim E_J^2/E_C$, hence the visible reduction shown in Fig. (a). This is important, as we can now argue from a graphical standpoint why the charge dispersion increases for higher excited states: the splitting size decreases and $$ \frac{d}{d n_g} \bigg((n-n_g)^2 \bigg) $$ increases with respect to the fixed values of $n$ that define the relevent energy level. 66