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$.