Extrapolation of solvability of the parabolic $L^p$ Neumann problem on bounded Lipschitz cylinders

A recent result of the first author with Li and Pipher has established the extrapolation of solvability of the $L^p$ parabolic Neumann problem on unbounded graph domains of the form $\Omega=\{(x',x_n):\,x_n>\varphi(x')\}\times\mathbb R$, where $\varphi:\mathbb R^{n-1}\to\mathbb R$ is a Lipschitz function. The result shows that under the assumptions that the $L^p$ parabolic Neumann problem for the equation $Lu=-\partial_t u+\mbox{div}(A\nabla u)=0$ in $\Omega$ and also the $L^{p'}$ parabolic Dirichlet problem for the adjoint equation $L^*u=\partial_t u+\mbox{div}(A\nabla u)=0$ in $\Omega$ are solvable, then also the $L^q$ parabolic Neumann problem for the equation $Lu=0$ in $\Omega$ is solvable for all $1<q<p$.
However the mentioned paper does not answer the question whether the same claim is also true for domains of the form $\mathcal O\times\mathbb R$, where $\mathcal O$ is a bounded Lipschitz domain (in spatial variables) since this case does not follow from our argument for the unbounded case. Indeed, the bounded Lipschitz cylinder case requires a significantly different approach which we present in this article and establish an analogous result when $\mathcal O$ is a bounded Lipschitz domain.