Well, lots of time have passed and now I have an explicit formula for numerosity. It gives the same (up to an infinitesimal)  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 the 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}}\simeq \frac12$ 

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 {3 k^2 + k } 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)\simeq \frac1{\sqrt{2}}$

This is infinitesimally close to your result.

Perhaps before calling a surreal number a numerosity we should exclude the ifinitesimal part from it, which in surreal numbers can be done via decomposing the number into Conway normal form. If we do so, our results will completely coincide.

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