Parabolic Regularity with Neumann B.C Consider the parabolic problem in the cylinder of base $B$, the unit ball,
$$
\partial_t u -\text{div}\left( A(x) D u +F(t,x)\right)=0 \text{ in } (0,T)\times B,
$$
with $(ADu +F)\cdot \nu=0$ on $(0,T)\times \partial B$, and $u(t=0)=0$,in 
Suppose $A$ is smooth (let us say $C^{2}$), symmetric, $2 I_d \geq A\geq I_d$.    What is the regularity required on $F$ for $u$ to be bounded independently of $T$? I am looking for an estimate of the form 
$$ 
\sup_{(0,T)\times B} |u(t,x)| \leq C \| F \|_{L^p((0,T);H^{q}(\text{div},B))} 
$$ 
where $H^q(div,B)$ means the space of functions in $L^q$ with divergence in $L^q$ (thank you Hannes),  where $C$ depends on $A$ (potentially) but not on $u$, $F$, or $T$. 
 A: Let's see how this goes. First, let me say that the continuity estimate you are looking for is contained in [3, Theorem 4.5] under the integrabilities as I had already mentioned, so $q > n$ and $p > 2(1-n/q)^{-1}$. Unfortunately, to get the uniformity w.r.t. $T$ of the constants, we have to take a closer look at the argument and do a slight detour. If everything there is already clear for you, you might probably skip some of this..
I am going to leave out the "$(B)$" declarations for the function spaces.
Doing the usual integration by parts thing and recognizing that $A$ is bounded over $B$, we see that the weak form of your problem is given by $$u'(t)-\nabla\cdot A\nabla u(t) = \nabla \cdot F(t) \quad \text{in}~\bigl(W^{1,q'}\bigr)' \quad \text{f.a.a.}~t \in J$$ with $$\langle \nabla\cdot G,\varphi\rangle = \int_B (G,\nabla \varphi)_{\mathbb{R}^n}$$ and $J$ being the time interval, plus the initial value. Then $u$ is a solution to this equation if and only if $v(t) := e^{-t}u(t)$ is a solution to $$v'(t)-\nabla\cdot A\nabla v(t) + v(t) = \nabla \cdot e^{-t}F(t) \quad \text{in}~\bigl(W^{1,q'}\bigr)' \quad \text{f.a.a.}~t \in J\tag{1}$$
with $v(0) = 0$.
For this equation, we can appeal to maximal parabolic regularity of $-\nabla \cdot A \nabla + 1$ over $W^{-1,q} = (W^{1,q})'$ [2, Sect. 11]/[3, Thm. 4.6] with its domain $W^{1,q}$, which gives a unique solution $$v \in W^{1,p}(J;W^{-1,q}) \cap L^p(J;W^{1,q}) \hookrightarrow C(\overline J;(W^{-1,q},W^{1,q})_{1-1/p}) \tag{2}$$ together with the estimate $$\|v\|_{W^{1,p}(J;W^{-1,q}) \cap L^p(J;W^{1,q})} \leq C_S \|\nabla\cdot e^{-\cdot}F(\cdot)\|_{L^p(J;W^{-1,q})}.$$
(See [1, Thm. III.4.10.2] for the embedding and Chapters III.1.5 and III.4.10 for more. Note also that the reasoning also works for the abstract domain $\mathcal{A}_q$ instead of $W^{1,q}$ as in [3].) Here $C_S$ is the operator norm of the parabolic solution operator $(\partial -\nabla\cdot A \nabla + 1)^{-1}$. Now it remains to observe that [3, Lem. 4.8]
$$(W^{-1,q},W^{1,q})_{1-1/p} \hookrightarrow C(\overline B) \quad \text{if}~q>n~\text{and}~p>2(1-n/q)^{-1}, \tag{3}$$
so together with (2)
$$\|v\|_{C(\overline J \times \overline B)} \leq C_{IP}C_{MR}C_S\|\nabla\cdot e^{-\cdot}F(\cdot)\|_{L^p(J;W^{-1,q})},$$
where $C_{IP},C_{MR}$ and $C_S$ are the norms of the interpolation space embedding (3), of the maximal regularity space embedding (2) and of the solution operator. 
Going back to $u(t) = e^tv(t)$, we thus have $$\|u\|_{C(\overline J \times \overline B)} \leq C_{IP}C_{MR}C_S\|\nabla\cdot F\|_{L^p(J;W^{-1,q})},$$ as desired.

Now to the dependence of the collected constants on the interval. Clearly, $C_{IP}$ is independent of $|J|$. For the others, this is in general not true, but luckily the semigroup generated by $-\nabla\cdot A\nabla + 1$ admits a negative growth bound (here the "$+1$" becomes really useful), so not only finite intervals $J= [0,T]$ are admissible in the foregoing stuff, but also $J = \mathbb{R}^+$, see [4].
Given $0 < T < \infty$, prolong $\nabla \cdot e^{-\cdot} F(\cdot)$ by zero and solve (1) on $\mathbb{R}^+$ with the solution $v^\infty$, then $v^\infty_{\restriction [0,T]}$ is the solution $v^T$ to (1) on $[0,T]$ and
$$\begin{align}\|v^T\|_{W^{1,p}(0,T;W^{-1,q}) \cap L^p(0,T;W^{1,q})} & \leq \|v^\infty\|_{W^{1,p}(\mathbb{R}^+;W^{-1,q}) \cap L^p(\mathbb{R}^+;W^{1,q})} \\ & \leq C_S^\infty \|\nabla \cdot e^{-\cdot} F(\cdot)\|_{L^p(\mathbb{R}^+;W^{-1,q})} & \\ & = C_S^\infty \|\nabla \cdot e^{-\cdot} F(\cdot)\|_{L^p(0,T;W^{-1,q})},\end{align}$$
so $C_S^T \leq C_S^\infty$ since $C_S^T$ was an operator norm. 
For $C_{IP}$ as in (2), see Lemma III.4.10.1 (ii) and Remark III.1.4.3 in the book of Amann [1].

[1] Amann, Herbert, Linear and quasilinear parabolic problems. Vol. 1: Abstract linear theory, Monographs in Mathematics. 89. Basel: Birkhäuser. xxxv, 335 p. (1995). ZBL0819.35001.
[2] Auscher, Pascal; Badr, Nadine; Haller-Dintelmann, Robert; Rehberg, Joachim, The square root problem for second-order, divergence form operators with mixed boundary conditions on $L^p$, J. Evol. Equ. 15, No. 1, 165-208 (2015). ZBL1333.47034. 
[3] Disser, Karoline; Rehberg, Joachim; terElst, Tom, Hölder estimates for parabolic operators on domains with rough boundary, accepted in Ann. Sc. Norm. Sup. Pisa.
[4] Dore, Giovanni, $L^p$ regularity for abstract differential equations, Komatsu, Hikosaburo (ed.), Functional analysis and related topics, 1991. Proceedings of the international conference in memory of Professor Kôsaku Yosida held at RIMS, Kyoto University, Japan, July 29-Aug. 2, 1991. Berlin: Springer-Verlag. Lect. Notes Math. 1540, 25-38 (1993). ZBL0818.47044.
