  We discuss the problem of classifying polynomials $p : \mathbb R^2_+
\rightarrow (0, \infty)$ for which $\frac{1}{p}=\{\frac{1}{p(m, n)}\}_{m, n
\geq 0}$ is joint completely monotone, where $p$ is a linear polynomial in $y.$
We show that if $p(x, y)=a+b x+c y+d xy$ with $a > 0$ and $b, c, d \geq 0,$
then $\frac{1}{p}$ is joint completely monotone if and only if $a d - b c \leq
0.$ We also present an application to the Cauchy dual subnormality problem for
toral $3$-isometric weighted $2$-shifts.
