Constant curvature difference surfaces Surfaces of Constant Mean Curvature CMC and constant Gauss Curvature Product $K$ are well-known in differential geometry of surfaces of revolution. Denote half sum and half difference curvatures as
$$ k_1+ k_2 = 2 H_s;  \, k_1- k_2 = 2 \, H_d. $$
Significance of $H_d$
Now I assume that  $H_d$ may as well be of equal interest with fundamental importance. The assumption is  motivated by its representation in Mohr's Circle radius of curvature for stress, shell curvatures, moment of inertia (in Figure) among other such tensors. As is known, material stress failure theories in structural mechanics operate on stress/strain and other tensoral differences as shear entities. It occurs as an important curvature invariant in the equation of Mohr's Circle :
$$\boxed{ (k_n-H_s)^2 +\tau_g^2= H_d^2} $$

I made a search in some textbooks that could be accessed. I was not lucky in finding references for the principal curvature difference profiles. Also the surface cannot be viewed as a particular case of CMC surfaces due to its sign change of curvature and so it is distinctly different. The  calculation has a different expression and plots differently.
If $\phi$ is slope of tangent to meridian, primes on meridian arc
Const $H=H_s$ Mean curvature Delaunay meridians
$$\phi^{'}+\frac{\cos \phi}{r} =2 H_s$$
With initial radius $r=r_1$ at $\phi=0$
$$ \cos \phi = \frac{H_s (r^2-r_1^2)+r_1}{r}$$
$$(-k_1,k_2)= (H_s(r^2+r_1^2)-r_1,\, H_s(r^2-r_1^2)+r_1)\,$$
Const $H_d$ Difference curvature meridians
First Order ODE:
$$-\phi^{'}+\frac{\cos \phi}{r} = 2 H_d \tag1$$
First order ODE has the disadvantage of numerically computing indefinitely below points of singularity on $r=0$ axis so there appear multiple spindle meridian profiles below this symmetry axis.
Second Order ODE:
$$\phi^{''}= 2 H_d \tan\phi\, (\phi^{'} +2 H_d) \tag2 $$
$$ \cos \phi = \frac{r}{r_1}+2\, H_d\, r\, log \,\frac{r}{r_1}\tag3 $$
$$(k_1,k_2)=(\frac{1}{r_1}+2H_d(1+log \,\frac{r}{r_1}),\frac{1}{r_1}+2H_d \,log \,\frac{r}{r_1} ) \tag4 $$
$$ @\,r=0,\phi \rightarrow \pi/2$$
Second order ODE has the advantage of stopping computation at point of singularity so there appears no profile below symmetry line $r= 0.$
Shown below are profiles of constant difference $H_d$ of principal curvatures.
Three distinct shapes occur.
Progressive loops $ H_d <-0.5$ above $r=1$ ;
Two types below $r=1$
Ovaloids between cylinder and sphere $ 0>H_d>-0.5;$ and,
Profiles with Inflection Point  for $ H_d >0 $ occurring at $ r=r_1e^{-(1+1/(2 r_1H_d))}. $
All profiles meet the axis of symmetry normally, however these are not umbilical points.
$ H_d$ Unduloids" />
Thanks in advance for your comments and for any references available on the topic.
 A: I see a somewhat different description of the surfaces of revolution with $H^2-K$ equal to a positive constant, say $H^2-K=\delta^2$, where $\delta>0$ is constant.  Note that scaling a surface by a constant $\lambda\not=0$ produces a new surface with $H^2-K$ divided by $\lambda^2$, so the actual magnitute of $\delta>0$ is not important for the possible shapes of such surfaces.  Thus, your pictures strike me as somewhat misleading.  
Instead, the key factor determining the shape of the corresponding surface of revolution is another constant, which can be described as follows:  First, let $\bigl(f(s),g(s)\bigr)$ be a unit speed curve in the plane, and consider the surface of revolution parametrized by
$$
X(s,\theta) = \bigl(f(s)\,\cos\theta,\ f(s)\,\sin\theta,\ g(s)\bigr).
$$
One easily computes that the principal curvatures of this surface are 
$$
\kappa_1 = g'(s)f''(s)-f'(s)g''(s)\quad\text{and}\quad \kappa_2 = -g'(s)/f(s),
$$
where, of course, $f'(s)^2+g'(s)^2 = 1$.  Then the equation we want to study is
$$
g'f''-f'g'' + g'/f = 2\delta.
$$
Using the relation $f'f''+g'g'' = 0$, this becomes the relations
$$
f'' = g'(2\delta-g'/f)\quad\text{and}\quad g'' = -f'(2\delta-g'/f).
$$
Scaling the surface by $\lambda$ induces the scaling
$$
(f,g,s,\delta) \mapsto (\lambda f,\lambda g, \lambda s, \delta/\lambda),
$$
so we can reduce to the case $\delta=\tfrac12$ without loss of generality.  Set $h=g'$
and we have the equation $h' = -f'(1-h/f)$ or $fh'-hf'+ff'=0$, which can be made exact by dividing by $f^2$, yielding the first integral
$$
\frac{h}{f} + \log|f| = C,
$$
so $h = f(C-\log|f|)$, and, since $(f')^2 = 1-h^2$, we must have $|f|(C-\log|f|) \le 1$.  In particular, we have
$$
\frac{(\mathrm{d}f)^2}{\bigl(1-f^2(C-\log|f|)^2\bigr)} = (\mathrm{d}s)^2,
$$
and, since $(\mathrm{d}g)^2 = h^2\,(\mathrm{d}s)^2$, we arrive at the relation
$$
\bigl(1-f^2(C-\log|f|)^2\bigr)\,(\mathrm{d}g)^2 - f^2(C-\log|f|)^2\,(\mathrm{d}f)^2 = 0
$$
It is the constant $C$ that determines the shape of the resulting profile curve and hence the surface.  Obviously, the null curves of this quadratic form have to lie inside the region where $f^2(C-\log|f|)^2 \le 1$, which are bounded by lines defined by
$|f|(C-\log|f|) = \pm1$.  
For example, when $C>1$, there are three functions $\rho_i$ on $(1,\infty)$ such that
$$
0 < \rho_1(C) < \mathrm{e}^{C-1} < \rho_2(C) < \mathrm{e}^C < \rho_3(C)
$$
such that the function $1-f^2(C-\log|f|)^2$ is nonnegative on the three intervals, 
$(-\rho_3(C),-\rho_2(C))$, $(-\rho_1(C),\rho_1(C))$, and $(\rho_2(C),\rho_3(C))$.  Consequently, there wil be two distinct solutions $(f,g)$:  One where $f$ and $g$ are periodic and $|f|$ is bounded by $\rho_1(C)$, and another where $f$ is periodic, oscillating between a minimum of $\rho_2(C)$ and a maximum of $\rho_3(C)$.
Meanwhile, when $C<1$, there is a single function of $C$, say $\rho_4(C)>\mathrm{e}^C$ such that the function $1-f^2(C-\log|f|)^2$ is nonnegative on the interval $(-\rho_4(C),\rho_4(C))$, and there will be a single solution (up to translation and reflection in $g$) for which $f$ lies in this interval.
There is a similar story for $C=1$, except that there is an exact solution, $(f,g)=(1,s)$ (the cylinder of radius $1$) and then two solutions, one with $0<f<1$ that asymptotically approaches the cylinder from the interior but pinches to a point on the axis of rotation in finite distance, and one with $1<f<M$ (where $M\approx3.591121477$ satisfies $M(1-\log M) = -1$) that is asymptotic to the cylinder of radius $1$ in both directions, but bulges out to a maximum diameter of $M$ at its 'equator'.  The resulting surface of revolution has a circle of self-intersections. 
It is probably worth noting that, near $f=0$, i.e., the axis of rotation, the function $g$ is not a $C^2$ function of $f$ (obvious, since, otherwise, the surface would have an umbilic there), but has an expansion of the form
$$
g = g_0 \pm \tfrac14 f^2\bigl(2C+1-\log|f|\bigr) + O(f^3)
$$
Thus, the surface is only $C^1$ when it intersects the axis of rotation. 
Note that all of these solutions have $\delta=\tfrac12$, but they have quite different behaviors.
A: ..a small addition to Professor Bryant's answer
For this class of surfaces, where $k1-k2 = constant$, we have introduced the name "Costant Skew Curvature Surfaces" (CSkC-surfaces) and we have studied an aspect concerning the Bonnet-surfaces.
It is well know that the famous question that Bonnet asked was: "When does there exist an isometric embedding $x:M^2 \rightarrow R^3$ such that the mean curvature function of the immersion is $H$?" 
Our work was born from the question: Can a surface be CSkC and Bonnet at the same time, and, if that is the case, what does it represent?
We showed that the CSkC-surfaces with principal curvatures ($k_1$ and $k_2$) nonconstant cannot contain any Bonnet-surfaces, so if and only if $k_1$ and $k_2$ are both constant the class of CSkC-surfaces can admit Bonnet-surfaces.
This means that the only CSkC-surfaces for which exists a nontrivial isometric deformation preserving the mean curvature $H$ are (patch of) a circular cylinders.
see (Alexander Pigazzini, Magdalena Toda): https://projecteuclid.org/euclid.jgsp/1518577293
Another aspect of the CSkC-surfaces is that if we setting $k_1-k_2=constant$ renders the shape equation for an elastic membrane equivalent to the Schrodinger equation for a particle on the same surface.. then the same equations have the same solutions...
see (Victor Atanasov, Rossen Dandoloff):   https://iopscience.iop.org/article/10.1088/0143-0807/38/1/015405/pdf
By relating the two works we can for example say that:
"if you want to look for an isometric deformation that preserve this relation (shape equation equivalent to the Schrodinger equation) then you can get the result only if the surface is (patch of) circular cylinder, and then search through its isometric deformations".
