Well, lots of time have passed and now I have an explicit formula for numerosity. In case of uniform lattices it gives the same differences between numerosities of sets as your formula, but can express the numerosities in a precise form.
There is a clear way to express numerosities of sequences via $\omega$, which is the germ of the identity function at infinity or half the numerosity of $\mathbb{Z}$. It also can be considered a surreal number, given the canonical embedding of Hardy fields into surreals. This measure is more accurate than asymptotic density, because it gives the exact value.
Suppose you have a strictly increasing sequence $a_k\ge0$, where $k\in\mathbb{Z}, k\ge0$.
To find the numerosity, you have to apply to your sequence the operator $N(a_k)=\left(D\Delta^{-1}a_k\right)^{[-1]}(\omega)$, where $f^{[-1]}$ is the inverse function.
The following Wolfram Language code does the thing:
a[k] := k^2
SolveValues[D[Sum[a[k], k], k] == \[Omega], k] /. C[1] -> 0 //
FullSimplify // Expand
Inverse code to find a sequence with desired numerosity:
S = Log[\[Omega]];
DifferenceDelta[Integrate[Normal[SolveValues[S == k, \[Omega]]], k],
k] /. C[1] -> 0 // Last // FullSimplify // Expand
Basically, we represent a surreal number as a germ at infinity, the germ as a divergent integral, then divide the integral into segments of area $1$, and the centers of mass of these segments are the set of the desired numerosity.
Let us apply the formula to examples from your question and see what comes.
- $N(\{0,2,4,6,...\})-N(\{1,3,5,7,...\})=N(2k)-N(2k+1)=\frac{\omega }{2}+\frac{1}{2}-\frac{\omega }{2}=1/2$
This coincides with your result.
- $N(\{0,1,4,9,16,...\})-N(\{0,2,6,12,20,...\})=N(k^2)-N(k(k+1))=\sqrt{\omega +\frac1{12}}+\frac{1}{2}-\sqrt{\omega +\frac{1}{3}}$
This is infinitesimally close to your result.
- $N(\{0,1,5,12,22,...\})-N(\{0,2,7,15,26,...\})=N({\frac {3k^{2}-k}{2}})-N(\frac {k (3 k + 1) } 2)=\frac{1}{3} \sqrt{6 \omega +1}+\frac{2}{3}-\left(\frac{1}{3} \sqrt{6 \omega +1}+\frac{1}{3}\right)=\frac13$
This coincides with your result.
- $N(\{0,1,3,6,10,..,\}) - \sqrt{2} N(\{0,1,4,9,16,...\})=N(\frac12k(k+1))-\sqrt{2}N(k^2)=\sqrt{2 \omega +\frac{1}{3}}-\left(\sqrt{2 \omega +\frac{1}{6}}+\frac{1}{\sqrt{2}}\right)$
This is infinitesimally close to your result.
Some other examples:
$N(1/3+k+k^2)=\sqrt{\omega}$
$N(k^4)=\frac{1}{30} \sqrt{30 \sqrt{900 \omega +30}+225}+\frac{1}{2}$
$N(7^k)=\log_7 \left(\frac{6 \omega }{\ln (7)}\right)$