I've managed to essentially brute force the problem by calculating the Hessian of the function, and showing that its determinant and trace are non-negative.
This was done by using a change of variable to reduce the problem to showing that two certain polynomials are positive over a subset of $[0,1]^2$, proving that it's non-negative in a neighborhood of its zeros, and numerically checking that it's positive away from them.
This solution feels a bit too messy for me, so I was wondering if there isn't a cleaner approach one could use. (I'm aware we could use Sylvester's criterion to simplify the numerical step, but I'd like to avoid using that as well if possible.)
For reference, the expression of the Hessian is. $$H(x,y) = \begin{bmatrix} -s(1-s)(1-2s)(s-t) + s^2(1-s)^2 && -s(1-s)t(1-t) \\ -s(1-s)t(1-t) && t(1-t)(1-2t)(s-t) + t^2(1-t)^2 \end{bmatrix}.$$
where $s=\operatorname{sigmoid}(x),\ t=\operatorname{sigmoid}(y)$.