Kernel bounds for parabolic operators having first-order degeneracy at the boundary

We study kernel estimates for parabolic problems governed by singular elliptic operators \begin{equation*} \sum_{i,j=1}^{N+1}q_{ij}D_{ij}+c\frac{D_y}{y},\qquad c+1>0, \end{equation*} in the half-space ${\mathbb{R}}_{+}^{N+1}=\{(x,y):x\in {\mathbb{R}}^N,y>0\}$ under Neumann boundary conditions at $y=0$.