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)$


  [1]: https://en.wikipedia.org/wiki/Indefinite_sum